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CHAPTER  1 
INTRODUCTION 


In  this  report  we  develop  a  general  statistical  model  for  the  analysis  of 
explosives  performance  data.  The  data  of  interest  consist  of  observations  of  an 
arbitrary  measure  of  performance  that  can  be  derived  from  the  transient 
pressure-induced  responses  of  piezo-electric  gages  located  in  the  vicinities  of 
underwater  explosions.  Applications  of  the  model  concern  the  variations  of 
performance  within  broadly  defined  classes  of  explosive  charges  and  the 
predictions,  and  comparisons,  at  specified  ranges,  of  the  performances  of  charges 
belonging  to  these  classes.  The  observations  are  complicated  by  the  presence  of 
gage  calibration  errors  and  various  errors  that  arise  from  the  measurement  and 
data  processing  techniques  employed. 

The  traditional  and  current  method  for  analyzing  data  of  this  kind  (see 
Cole,  Reference  1,  p.  240)  is  to  express  the  relationship  between  the  scaled 

performance  variable  and  the  scaled  distance  from  the  charge  as  a  power  law 

\ 

whose  parameters  are  determined  from  an  ordinary  least  squares  fit  of  a 
straight  line  to  the  logarithmically  transformed  data.  Relationships  obtained 
in  this  manner  are  popularly  called  "similitude  equations"  from  the  principle  of 
similarity  that  underlies  the  scaling  of  explosion  shock  wave  phenomena.  The 


*Cole,  R.  H. ,  Underwater  Explosions  (New  Jersey:  Princeton  Univ.  Press,  1948). 
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usual  method  of  scaling  is  called  Hopkinson  or  cube-root  scaling  (see  e.g. 

Snay,  Reference  2),  which  refers  to  the  fact  that  the  time  and  length  scale 
factors  are  proportional  to  the  cube  root  of  the  charge  weight  (or 
any  other  proportional  measure  of  the  explosion  energy).  Recently  Goertner 
(Reference  3)  has  extended  this  scaling  method  to  include  variations  with 
ambient  water  sound  speeds  and  densities.  It  is  now  customary  to  compile 
similitude  equations  for  a  large  variety  of  explosive  classes  for  the  following 
measures  of  performance:  the  peak  pressure,  a  characteristic  time  constant,  and 
the  impulses  and  energies  per  unit  area  delivered  to  a  given  location  by  the 
shock  wave  within  various  multiples  of  the  time  constant. 

The  application  of  more  sophisticated  statistical  techniques  to  this  area 
has  been  blocked,  perhaps  primarily,  by  complicated  dependencies  within  the 
data.  Correlations  exist  among  the  observations  obtained  from  a  single  test,  as 
pointed  out  by  Brown  (Reference  4),  and  among  observations  obtained  with  the 
same  gages  and  gage  calibration  constants.  The  mixed  linear  models  necessary 
for  an  adequate  statistical  treatment  of  such  data  constitute  an  active  area  of 
current  research. 

The  model  that  we  develop  below  generalizes  the  current  approach  by 

/ 

accepting  arbitrary  measures  of  performance  and  allowing  regression  functions  of 
arbitrary  form  that  are  linear  in  their  coefficients  (such  as  higher  degree 

^Snay,  H.  G. ,  "Model  Testsand  Scaling,"  ROLTR  63-257,  1  Dec  1964. 

^Goertner,  J.  F. ,  "Scaling  Underwater  Explosion  Shock  Waves  for  Differences  in 
Ambient  Sound  Speed  and  Density,"  NSWC  TR  80-491,  18  Dec  1980. 

^Brown,  R.  H. ,  "Analysis  of  Data  When  Several  Sources  of  Variation  are 
Present,"  Explosives  Research  Memorandum  22,  Navy  Dept.  Bureau  of  Ordnance, 

1  Dec  1944. 


NSWC  TR  82-74 


polynomials).  As  ch«  principle  of  shock  wave  similarity  is  still  adhered  to  we 
continue  to  refer  to  the  equation  of  the  regression  mean  is  the  similitude 
equation  for  the  particular  explosive  class— a  term  that  we  will  use  to 
emphasize  the  fact  that  explosive  charges  are  more  properly  viewed  as  members  of 
a  class  of  objects  that  in  many  respects  are  alike  but  which  differ  in  ways  that 
affect  the  observed  measures  of  performance. 

The  model  extends  the  presently  used  approach  by  explicitly  including 
sources  of  random  variation  in  its  formulation.  In  the  interest  of  model 
simplicity  our  philosophy  has  been  to  include  only  those  sources  that  are 
thought  to  produce  significantly  large  effects  and  are  unavoidable. 

Complicating  effects  that  are  avoidable  or  correctible  will  be  presumed  to  have 
been  eliminated  either  by  an  appropriate  reprocessing  of  the  data  or  by 
modifications  of  the  experimental  techniques. 

A  possible  model  deficiency  is  that  no  explicit  treatment  of  so  called 
batch  effects  is  included.  This  refers  to  well  recognized  performance 
variations  among  charges  taken  from  different  batches  or  preparations  of  the 
same  explosive  material.  It  was  felt  that  batch  effects  did  not  justify  the 
further  complication  of  an  already  complex  model  and  that  they  could  be  handled 
in  another  manner  such  as  by  treating  different  batches  as  different  explosive 
classes,  by  increasing  the  number  of  batches  and  randomizing  the  charge 
selection,  and  by  improving  the  explosive  preparation  quality  control.  In  any 
case  the  use  of  the  model  should  be  made  with  the  possibility  of  batch  effects 
borne  in  mind,  and  a  thorough  examination  of  the  model  residuals  for  the 
presence  of  these  and  any  other  systematic  effects  is  recommended  (see  Section 
5-1). 


-3 


The  report  is  divided  into  five  chapters.  Chapter  2  discusses  the 


development  of  the  model  in  detail.  The  model  may  be  described  as  an 
intra-class  regression  model  with  an  additive  error  term  consisting  of  a  gage 
class  calibration  error,  a  within-class  performance  variation,  and  a  general 
experimental  error.  Under  a  suitable  transformation  of  the  performance 
variable,  multivariate  normality  of  the  errors  is  assumed.  Chapter  3  discusses 
the  estimation  of  the  model  parameters  by  the  methods  of  maximum  likelihood  and 
restricted  maximum  likelihood.  The  derivation  of  derivatives  needed  for  an 
iterative  solution  of  the  likelihood  equations  follows  the  approach  of  Harville 
(Reference  5).  In  Chapter  4  we  find  a  brief  description  of  both  unconstrained 
and  constrained  Newton-Raphson  and  method  of  scoring  optimization  techniques  and 
related  topics.  And  finally  in  Chapter  5  we  give  some  practical  applications  of 
the  model. 


^Harville,  D.  A.,  ..a  imi  Likelihood  Approaches  to  Variance  Component 

Estimation  and  to  K_j.ated  Problems,"  Journal  of  the  American  Statistical 
Assoc.,  Vol.  72,  No  358,  1977,  p.  320. 
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CHAPTER  2 


MODEL  DEVELOPMENT 


We  will  denote  the  response  and  regressor  variables  of  the  model  as  y  and  x 
respectively,  often  with  subscripts  to  specify  a  particular  observation. 

Vectors  in  the  model  will  be  indicated  by  use  of  the  underbar  notation;  hence,  a 
sample  of  response  variables  will  appear  as  £.  No  notational  distinction  will 
be  made  between  realized  and  random  samples,  but  this  difference  should  be 
apparent  from  the  context. 


In  the  theory  of  linear  models  it  is  common  to  deal  with  transformed 
response  and  regressor  variables  to  promote  variance  homogeneity,  model 
simplicity,  and  other  desirable  model  properties.  Thus,  we  define  y  as  a 
possibly  transformed  value  of  a  scaled  measure  of  performance,  as  discussed 
earlier,  and  x  as  a  possibly  transformed  scaled  distance  taken  to  be  zero  at  the 
charge  center.  This  is  in  accord  with  past  derivations  of  shock  wave  similitude 
equations  in  which  straight  lines  are  fitted  to  the  logarithms  of  the  scaled 
data.  Hence,  the  present  model  will  be  compatible  with  these  forms. 


Development  of  the  model  will  be  based  upon  a  number  of  reasonable 
assumptions  which  will  appear  throughout  this  section.  To  these  and  the 
assumption  of  shock  wave  similarity  already  made  we  add  that  the  water  between 
the  charge  and  gages  is  assumed  t'  be  homogeneous  so  that  disturbances  are 
propaged  through  the  water  in  a  regular  manner,  and  we  assume  the  values  of  x 
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to  be  accurately  known.  Recent  tests  of  the  last  assumption  have  verified  its 
accuracy  when  distances  are  determined  from  the  measured  times  of  arrival 
(Reference  6). 


2-1  BASIC  MODEL  FOR  A  SINGLE  OBSERVATION 


In  the  sample  of  y  values  let  y £ jjcm  be  Che  observation  made  in  the  jth 
shot  of  the  ith  explosive  class  with  a  gage  of  the  kth  gage  class  having 
calibration  index  m.  The  gage  is  located  at  the  (transformed,  scaled)  distance 
Xijkm*  T*ie  8a®e  class  index  refers  to  one  of  several  broad  classes  of  gages 
such  as  1/4  inch  tourmaline,  3/8  inch  tourmaline,  1/2  inch  tourmaline  etc.  In  a 
typical  test  it  is  common  to  employ  a  string  of  perhaps  10  to  12  gages  placed  so 
that  smaller  diameter  gages  are  grouped  closer  to  the  charge  and  larger  diameter 
gages  grouped  farther  from  the  charge.  Typically  gages  from  3  to  4  gage  classes 
are  used.  A  gage  is  usually  recalibrated  prior  to  each  test  program  and 
assigned  a  gage  constant  (units  of  picocoulombs/psi)  which  is  used  to  calculate 
all  pressures  measured  by  the  gage  until  it  is  recalibrated.  From  this  it  is 
clear  that  measures  of  performance  derived  from  the  various  pressure-time 
records  of  a  particular  gage  could  be  correlated,  i.e.,  all  affected  in  the  same 
manner.  During  the  lifetime  of  a  gage,  lasting  only  a  single  shot  to  perhaps 
several  years,  it  might  be  recalibrated  as  many  as  3  to  10  times. 


The  system  by  which  observations  will  be  indexed  is  described  as  follows. 
Let  C  and  K  be  the  total  numbers  of  explosive  classes  and  gage  classes  of 

^Gaspin,  J.  B.,  "Validation  of  a  Gage  Location  Method  for  Underwater  Explosion 
Tests,"  NSWC  TR  in  preparation. 
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interest,  respectively,  and  let  J.  be  the  total  number  of  shots  of  the  ith 
explosive  class.  Then  index  i  will  simply  run  from  1  to  C,  index  j  will  run 
from  1  to  J.,  and  index  k  will  run  from  1  to  K.  The  numbering  assigned  to 
particular  classes  or  shots  is  arbitrary.  In  the  same  sense  the  gage 
calibration  index  m  is  also  arbitrarily  assigned,  but  it  will  require  a  somewhat 
lengthier  explanation.  Basically,  we  wish  to  identify  those  observations 
obtained  from  the  same  gage  and  one  of  its  particular  gage  calibration  constants 
and  assign  to  these  observations  the  same  value  of  m.  We  will  establish 
separate  sets  of  m  values  for  the  observations  associated  with  each  gage  class 
of  interest.  For  the  kth  gage  class,  the  total  number  of  such  m  values  is 
equal  to  the  total  number  of  unique  ordered  pairs  (n,d)  among  the  observations 
associated  with  the  kth  gage  class,  where  n  is  a  unique  identifying  number  of 
the  gage  and  d  is  the  date  on  which  the  calibration  session  was  conducted, 
hence,  for  the  kth  gage  class  the  calibration  index  m  will  run  from  1  to  M^. 

Note  that,  although  several  measurements  can  have  the  same  indices  k  and  m,  an 
observation  is  uniquely  labeled  by  the  set  of  indices  (i,j,k,m),  since  only  a 
single  observation  can  be  obtained  from  a  particular  gage  on  a  particular  shot. 

Consider,  now,  the  jth  test  of  the  ith  explosive  class.  We  assume  that  the 
passage  of  the  shock  wave  through  the  water  induces  a  particular  unknown 
functional  relationship  between  the  quantity  we  are  seeking  to  measure  and  the 
travel  distance  x.  We  will  denote  a  representation  of  this  function  as  fj.j(x) 
for  the  jth  test  of  the  ith  explosive  class.  We  can  write  the  observation 
obtained  with  a  gage  having  indices  k  and  m  as 

yijkm  “  fij(xijkm)  +  eijkm»  (2-1) 
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where  e. is  the  deviation  of  the  measured  value  from  the  true  value  as 
i  jKm 

indicated  in  Figure  1.  Values  of  y  measured  by  other  gages  and  at  other  values 
of  x  will  be  scattered  about  f..(x)  in  some  manner  that  we  will  now  consider. 


Figure  1.  Response  y  versus  distance  x  for  jth  shot 

We  will  assume  the  error  e.  to  be  a  random  variable  with  a  mean  of 

l  jxm 

zero.  This  requires  that  any  bias  in  e^^  introduced,  for  example,  by  an 
incorrect  gage  size  correction  or  other  systematic  measurement  effect,  be 
removed  prior  to  the  statistical  analysis.  Current  thinking  by  experimentalists 
is  that  this  is  not  an  unreasonable  assumption.  Furthermore,  there  are  two 
compelling  reasons  for  it  from  theoretical  grounds.  First,  it  keeps  the  model 
from  becoming  unmanageably  complex.  Second,  and  most  importantly,  the  presence 
of  bias  terms  would  make  the  usual  similitude  equations  unestimable.  That  is, 
they  would  have  no  unique  solutions.  An  explanation  for  this  is  given  below  in 
Section  3-1.  Thus,  in  the  model  we  require 
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It  is  convenient  and  consistent  with  past  practices  to  represent  the 
functional  relationship  f^j(x)  by  a  polynomial  in  x.  This  can  be  written  in  a 
somewhat  more  general  fashion  for  the  jth  shot  of  the  ith  explosive  class  as 

f  (x  )  -  ♦'  u+  ,  (2-5) 

ij  ljkm  ljkm  ij 

where  *^jjan  i-8  a  vector  function  of  the  regressor1  variable  x^^  (prime 

denotes  the  transpose),  and  ji+  is  a  vector  of  unknown  parameter  values  that 

is  conformable  with  In  the  traditional  case  of  similitude  equations  we 

have  ♦ljkm  *  (1»  xijkm^*  For  the  case  of  a  (p-l)th  degree  polynomial  we  have, 

of  course,  ♦!  ..  *  (1,  x.  ,  x? ..  ).  Generally,  for  underwater 

— ljkm  ljKxn  ljkm  i  j  km 

explosions  a  low  degree  polynomial  with  p  *  2  or  3  will  be  adequate.  For  other 

applications  where  the  choice  of  p  >  6  may  seem  more  appropriate  it  may  be 

preferable  to  define  £  in  terms  of  Chebyshev  or  orthogonal  polynomials  to 

ijkm 

avoid  problems  of  ill  conditioning  (see  Seber,  Reference  7,  p.  214).  The 
generality  of  expression  (2~5)  should  be  fully  appreciated.  In  addition  to 
functions  of  a  single  dimension  it  will  admit  the  use  of  multidimensional 
functions  (that  are  linear  in  their  coefficients)  as  well.  In  the  subsequent 
development  we  will  assume  i.£jjcn  to  consist  of  p  elements  which  could  be  any 
of  these  types. 

The  relationship  between  y  and  x  will  vary  from  one  explosive  class  to 
another  and  also  between  shots  within  a  particular  explosive  class.  Within- 
class  variations  are  considered  to  be  caused  by  charge  fabrication  and 
preparation  methods  that  are  difficult  or  impossible  to  control  precisely,  and 

^Seber,  G.  A.  F. ,  Linear  Regression  Analysis  (New  York:  John  Wiley  &  Sons, 
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also  by  naturally  occurring  chemical  and  physical  inhomogeneities  randomly 
located  within  the  explosive  mixtures.  Mathematical  models  for  these 
relationships  must,  consequently,  be  of  a  random  nature.  We  can  view  the 
relationships  between  y  and  x  generated  by  the  shots  of  a  particular  explosive 
class  as  a  family  of  curves  each  of  which  is  similar  to  the  one  illustrated  in 
Figure  1.  This  idea  is  illustrated  in  Figure  3. 


Figure  3.  Response  versus  distance  curves  for  an  explosive  class 


A  mathematical  model  for  the  behavior  of  an  explosive  class  can  be 
constructed  by  allowing  the  parameters  jitj  of  equation  (2-5)  to  be  random 
variables.  Mean  values,  variances,  and  covariances  of  the  parameters  will  then 
be  constants  of  the  explosive  class.  If  we  denote  the  mean  of  as 
we  can  then  write 


(2-6 


which  decomposes  jjTj  into  the  sum  of  a  vector  of  fixed  effects  ^  and 
a  vector  of  random  effects  with  zero  means  Inserting  equations  (2-6) 

into  (2-5)  we  obtain 

2-7 
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f  (X  )  ■  y  +  ♦'  a*  .  (2-7) 

lj  ijkm  ijkm  i  ijkm  *ij 

Note  that  f^j(x)  is  now  regarded  as  a  random  function  consisting  of  a  mean 
curve  and  a  random  variation  about  that  mean. 

By  inserting  (2-7)  into  equation  (2-4)  we  obtain  the  basic  model  for  the 
analysis  of  explosion  test  data.  It  is  convenient  to  write  this  as 

y. .  ”  ±! .  j i. 

ijkm  ijkm  l 

+  a*  +  *'  8*  *  »  (2-8) 

km  ijkm  — lj  ijkm 

where  the  second  line  is  the  random  "error"  term  with  zero  mean.  The  basic 
model  is  recognized  as  a  linear- regression  model  with  a  complicated  error 
structure.  It  may  also  be  referred  to  as  a  mixed  model  or  mixed  regression 
mode 1 . 

Equation  (2-8)  expresses  the  contention  that  a  single  observation  may  be 
viewed  as  a  realization  of  a  random  variation  about  the  mean  response  of  the  ith 
explosive  class.  Furthermore,  it  holds  that  the  random  part  consists  of  the  sum 
of  (1)  a  random  gage  calibration  term,  (2)  the  random  performance  variations  of 
("identical")  charges  within  the  explosive  class,  and  (3)  a  random  experimental 
error  term.  And  finally,  it  indicates  that  the  mean  (or  expected  value)  of  the 
observations  of  the  ith  explosive  class,  that  corresponds  to  what  is 
traditionally  called  the  explosive  similitude  equation  (expressed  in  terms  of 
the  transformed  variable  y)  is  given  by 


2-8 
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E(y  )  *  *'  u 
ijkm  *"i  jkm  1 


at  distance  x.  . 

ijkm 


2-2  THE  MODEL  FOR  A  VECTOR  OF  OBSERVATIONS 


(2-9) 


We  now  express  the  model  in  matrix  form  in  three  steps  by  first  writing 
the  model  for  the  measurements  of  a  particular  shot  (i.e.,  for  fixed  values  of 
i  and  j),  then  combining  these  to  form  the  model  for  the  ith  explosive  class, 
and  finally  combining  several  single  explosive  class  models  to  obtain  a  model 
for  multiple  explosive  classes. 

To  reiterate  our  previous  discussion  concerning  indices  we  are  interested 
in  a  total  of  C  explosive  classes  and  K  gage  classes.  Also  we  let  the  number 
of  shots  in  the  ith  explosive  class  be  and  the  number  of  gage  calibrations 
in  the  kth  gage  class  be  M^.  Hence,  the  indices  take  on  the  values  i  ■ 

1 ,2 , . . . ,C;  j  m  1,2,... , J^;  fc  *  1,2, ...,K^  and  m  *  1,2,... ,M^ .  In  addx txon 
we  will  denote  the  number  of  observations  in  the  jth  shot  of  the  ith  explosive 
class  by  n. .  and  the  total  number  of  observations  as  N. 


For  fixed  values  of  i  and  j,  let  y. .  be  the  vector  of  n. .  observations 
for  the  jth  shot  of  the  ith  explosive  class.  The  model  for  this  vector  of 
observations  can  be  expressed  as 


1.  .  -X.  u  ♦  U. .a*  +  X. ,fl*.  +  e*., 
ij  iJ  i  ij  ij 


(2-10) 


where  X. .  is  the  matrix  of  regressor  variables  (whose  rows  consist  of  the  s' 
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row  vectors)  and  ejj  is  the  vector  of  values  corresponding  to  £ij*  Here 

the  vector  a*  denotes  the  complete  set  of  a*^  values  of  interest  in  the  full 
model,  ordered  according  to  gage  class.  That  is 

a  “  (®11 » • • • » °1M , »®2l * *  * • * ®2M  » • • • »®K1 » • • • »®KM  ^ ' • 

12  R 

t 

The  matrix  Ihj,  however,  is  specific  to  the  jth  shot  of  the  ith  explosive 
class  and  consists  of  0's  and  l's.  Since  only  a  single  gage  is  associated  with 
a  particular  observation,  Ihj  will  have  only  one  1  per  row.  An  example 
showing  the  elements  of  equation  (2-10)  explicitly  appears  in  Figure  4. 


To  assemble  the  model  for  the  complete  set  of  observations  of  the  ith 

explosive  class  we  define  the  observation  vector  of  interest  as  £  »  (yj i  , 

Zi2> • • • »ZiJ. '  »  the  corresponding  error  vector  as  e^  ■  (ej.1  »£j.2» •  •  **ij. )  » 

and  the  corresponding  vector  of  random  performance  effects  as  bj  »  (£J{, 

Ji2»  •  •  •  ifiij. )  •  Furthermore  we  define  the  following  matrices:  Xi  •  (X£i, 

Xi2»  •  •  •  ) ,  Oi  *  (^il  i ^£2 » •  •  • » )  ,  and  Wj  *  Block  Diag  (Xj.1  ,^i2 » •  •  •  »^ij . ) 

li  l 

With  these  the  model  for  ^  can  be  written  as 


Xu  +  U  a*  +  W  b*  +  e*. 


(2-1 


i  i 


l  l 


Finally,  we  form  the  full  matrix  model  by  defining  in  analogy  with  the 

above  x  -  <Xl »Z2» • • • >ZC>' »  £*  *  <£*’  »£*’ , . . .,e*’) U«  (U},U£, . . . ,U^) ’ ,  and  W 
Block  Diag  (Wj ,W2 , . . .  jWq).  Also  we  introduce  ^  . ^)’and  X  *  Block 

Diag  (X, jXj , . . . ,XC) .  In  terms  of  these  quantities  the  full  matrix  model  can 


be  represented  ar 
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K' 


£  *  Xji_  +  Ua*  +  Wb*  +  e*.  (2-12) 

As  before,  an  example  showing  the  full  matrix  model  in  more  explicit  form  is 
found  in  Figure  5. 


r: 


Equation  (2.12)  is  the  matrix  counterpart  of  equation  (2.8).  The  notations 
employed  to  describe  the  various  kinds  of  terms  are  similar  in  both  cases.  In 
words,  (2.12)  states  that  the  sample  of  observations  £  can  be  viewed  as  a 
realization  of  a  random  vector  in  an  N  dimensional  sample  space  that  deviates 
from  the  mean  point  Xjj^  by  a  vector  sum  composed  of  a  vector  of  gage 
calibration  errors  Ua*,  a  vector  of  performance  variations  Wb*,  and  a  general 
experimental  error  vector  £*.  Here,  of  course,  X  and  W  depend  upon  the 
transformed  reduced  distances  associated  with  the  sample  of  observations  and  U 
is  a  matrix  of  O's  and  l's  that  links  the  observations  with  the  random  gage 
effects  a*. 


> 


Viewing  as  a  variable,  X^  defines  a  hyperplane  in  the  N  dimensional 
sample  space  within  which  the  unknown  true  mean  is,  by  construction  of  the 
model,  postulated  to  lie.  In  Chapter  3-we  will  be  concerned  with  estimating  the 
position  of  the  true  mean  on  this  hyperplane  by  means  of  various  statistical 
criteria  of  choice.  For  example,  the  ordinary  least  squares  method  selects  the 
orthogonal  projection  of  £  on  the  hyperplane  as  the  estimated  value  of  Xj^. 

The  more  precise  maximum  likelihood  and  restricted  maximum  likelihood  methods  of 
estimation,  that  are  developed  in  this  report,  require  explicit  representations 
of  the  distributions  of  the  random  effects  in  the  model.  These  we  now  consider. 


,A 


2-12 


Example  showing  explicit  model  for  complete  observation  vector 
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2-3  DISTRIBUTIONAL  ASSUMPTIONS  AND  THE  DATA  VARIANCE  COVARIANCE  MATRIX 


ie  is  reasonable  Co  assume  Chat  Che  sample  of  observacions  £  will  have  a 
mulCivariaCe  normal  disCribuCion.  Hence,  we  wriCe 


Z  ~  N(Xu.,  V)  , 


(2-13 


where  X£  is  Che  mean  vecCor  of  Che  disCribuCion  and  V  is  Che  N  x  N  variance 
covariance  maCrix  whose  sCrucCure  we  will  consider  below.  As  noCed  earlier  Che 
adequacy  of  chis  assumption  may  require  a  transformation  of  Che  measuremenCs . 

For  Che  usual  quanCiCies  of  inCeresC  (peak  pressure,  decay  consCanC,  impulse  per 
unic  area,  and  energy  per  unic  area)  experience  suggescs  Che  logarichm  as  Che 
appropriaCe  Crans forma Cion.  Use  of  a  logarichmic  CransformaCion  also  supporCs 
Che  assumpCion  of  error  componenC  addiCiviCy  as  represented  in  equaCion  (2-3) 
since  calibration  effecCs  are  essentially  mulciplicative  on  the  pressure  (see 
Reference  1,  p.  183).  As  in  the  case  of  standard  regression  theory  it  will  be 
possible  Co  critically  examine  this  and  ocher  model  assumptions  through  an 
analysis  of  Che  model  residuals. 


Implicit  in  the  normality  assumpCion  for  £  is  the  assumpCion  that  all 
random  effects  of  Che  model  are  individually  normally  distributed.  Because 

•it  ft  ft 

a  ,  b  ,  and  £  are  affected  by  unrelated  error  sources,  they  are  taken  to 
be  mutually  independent.  We  will  assume  the  random  effects  to  have  the 
following  multivariate  normal  distributions, 


a*  ~  N(0, T  a2) 


(2-14 


*See  footnote  1  on  page  1-1. 
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(2-15) 


b*  ~  N(0,T  a2) 

e*  ~  N(0,I  a2)  (2-16) 


rfiifirc  I*  ®  Di3^  CY^,*.«,Y^,y^,...,Y\-g,...,  Yj^  ,  •  •  •  ,  y^  )  j 


M, 


*K 


Yk  =  a-2  Var(akm) ,  m  »  l,2,...,Mk, 


©i  *  a~2  Var(#h),  j  *  1,2, ...,Ji. 


(2-17) 


(2-18) 

(2-19) 


(2-20) 


.  .  .  .  .  2 
Hence,  r  is  the  matrix  of  gage  calibration  error  variances  divided  by  a 

and  T  is  the  matrix  of  the  variances  and  covariances  of  the  random  performance 
,  .  .  ,  2 

variation  parameters  divided  by  a  •  Correlations  between  the  calibration 
errors  of  gages  with  different  gage  class  and  calibration  indices  k  and  m  are 
assumed  to  be  zero,  as  are  the  performance  variation  parameters  of  different 
shots.  The  use  of  variance  covariance  ratios  is  for  later  mathematical 
convenience. 


As  mentioned  earlier,  relations  (2-17)  and  (2-18)  indicate  that  the  gage 
calibration  error  variances  are  taken  to  depend  on  the  gage  class  index  k  only 
and  not  on  other  factors  upon  which  they  might  reasonably  depend  such  as  reduced 
range,  charge  weight,  peak  pressure,  the  shock  wave  decay  constant,  or 
individual  details  of  the  gages.  This  is  a  simplifying  assumption  that  is 
thought  to  be  reasonably  accurate  and  of  practical  value. 


5 
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Using  these  distributional  assumptions  and  the  independence  of  a  ,  b  , 

* 

and  £  we  can  specify  the  data  variance  covariance  matrix  V.  From  (2-12)  and 

the  rule  for  forming  the  variance  covariance  matrix  of  a  linear  function  of 
* 

random  variables  we  obtain 


V  =  H'  a2 


(2-21) 


where 


H  =  I  +  uru  +  WTW 


(2-22) 


It  is  noted  from  the  block  diagonal  structures  of  W  and  T  that  WTW  is 
also  block  diagonal  with  blocks  of  size  n^x  n^j.  The  number 

of  observations  per  shot  is  usually  around  10.  The  relatively  simple  structure 
of  H  will  be  exploited  below  in  the  task  of  computing  functions  of  H  ^ . 


*If  £  “  Ax  with  E(x)  *  £  and  Var(x)  •  Z,  it  is  easily  shown  that  E(y)  *  Au_ 
and  Var(£)  *  E [ (^  —  E(^))(£  -  E(£))'l  *  AIA' . 


2-16 
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CHAPTER  3 

ESTIMATION  OF  THE  MODEL  PARAMETERS 


We  require  a  general  estimation  method  that  will  produce  efficient 
estimates  of  the  model  parameters,  and  provide  a  means  for  drawing  statistical 
inferences  concerning  the  constants  which  characterize  the  larger  populations 
from  which  the  sample  is  drawn.  Of  the  approaches  most  suited  to  mixed  models, 
the  maximum  likelihood  method  is  the  more  conceptually  and  theoretically 
attractive.*  In  the  past,  choice  of  this  method  has  been  often  avoided  because 
of  the  relatively  heavy  computational  burden  it  involves;  but  interest  has  been 
stimulated  in  recent  years  by  the  development  of  faster  computers,  more  rapid 
computational  techniques,  and  vcrious  theoretical  advances  (see  RejVruices  10 
through  16  and  the  survey  paper  by  Harville,  Reference  5).  Of  the  latter,  the 


*See  Searle  (Reference  8  and  Reference  9,  p.  458)  for  discussions  of  the 
suitability  of  different  approach  to  mixed  model  estimation. 

®Searle,  S.  R. ,  "Topics  in  Variance  Component  Estimation,"  Biometrics ,  Vol.  27, 
1971,  p.  1. 

^Searle,  S.  R.,  Linear  Models  (New  York:  John  Wiley  &  Sons,  Inc.,  1971). 

^Hartley,  H.  0.,  and  Rao,  J.  N.  K.,  "Maximum-Likelihood  Estimation  for  the 
Mixed  Analysis  of  Variance  Model,"  Biometrika ,  Vol.  54.  1  and  2,  1967,  p.  93. 

ll-Hemmerle,  W.  J.,  and  Hartley,  H.  0.,  "Computing  Maximum  Likelihood  Estimates 
for  the  Mixed  A.  0.  V.  Model  Using  the  W  Transformation,"  Technometrics , 

Vol.  15,  No.  4,  1973,  p.  819. 

l^Hemmerle,  W.  J.,  and  Lorens,  J.  A.,  "Improved  Algorithm  for  the  W-Transform 
in  Variance  Component  Estimation,"  Technometrics ,  Vol.  18,  No.  2,  1976, 
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concept  of  restricted  maximum  likelihood  (REML)  estimators  will  be  of  interest 

in  the  present  effort.  The  REML  method,  which  was  extended  to  the  general  mixed 

model  by  Patterson  and  Thompson  (Reference  17),  effectively  corrects  the  ML 

estimators  of  variance  components  for  losses  in  degrees  of  freedom  due  to 
estimation  of  the  fixed  effects.  In  this  report  we  will  develop  formulae  useful 
for  the  calculation  of  both  ML  and  REML  estimates.  In  this  effort  we  will  adapt 
and  follow  closely  the  results  obtained  by  Harville  (Reference  5). 


3-1  MAXIMUM  LIKELIHOOD  AND  OTHER  ESTIMATORS  OF  THE  SIMILITUDE  EQUATION 
PARAMETERS  u. 


Statistical  estimators  of  the  fixed  effects  in  linear  models  are  all 
closely  related.  Hence,  we  will  use  a  discussion  of  the  maximum  likelihood 
estimator  of  to  motivate  a  brief  discussion  of  other  estimators  of  the  fixed 
effects.  In  this  manner  some  of  the  advantages  of  ML  or  REML  estimators  over 
the  ordinary  least  squares  (OLS)  estimators,  currently  used  to  obtain  estimates 
of  explosive  similitude  equation  parameters,  can  be  highlighted. 


l^Jennrich,  R.  I. ,  and  Sampson,  P.  F.,  "Newton-Raphson  and  Related  Algorithms 
for  Maximum  Likelihood  Variance  Component  Estimation,"  Technometrics ,  Vol.  18, 
No.  1,  1976,  p.  11. 

l^Harville,  D.  A.,  "Some  Useful  Representations  for  Constrained  Mixed-Model 
Estimation,"  Journal  of  the  American  Statistical  Assoc.,  Vol.  74,  No.  365, 
1979,  p.  200. 

i^corbeil,  R.  R. ,  and  Searle,  S.  R. ,  "Restricted  Maximum  Likelihood  (REML) 

Estimation  of  Variance  Components  in  the  Mixed  Model,"  Technometrics ,  Vol.  18, 
No.  1,  1976,  p.  31. 

l^Corbeil,  r.  r. ,  and  Searle,  S.  R. ,  "A  Comparison  of  Variance  Component 
Estimators,"  Biometrics ,  Vol.  32,  1976,  p.  779. 

5see  footnote  5  on  page  1-4. 

^Patterson,  H.  D. ,  and  Thompson,  R.,  "Recovery  of  Inter-Block  Information 
When  Block  Sizes  are  Unequal,"  Biometrika ,  Vol.  58,  No.  3,  1971,  p.  545. 
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The  likelihood  function  L  is  defined  as  the  probability  of  the  observed 
sample  under  the  assumption  that  its  random  behavior  is  governed  by  some 
particular,  specified  family  of  distributions.  In  the  present  case  we  have 
assumed  that  £  has  a  multivariate  normal  distribution  with  unknown  mean  Xjj^  and 
variance  covariance  matrix  V.  Hence, 

L  =  (2ir)“N/2  |V|"1/2  exp[-(^-XVi),v'1(i-XvL)/2],  (3-1) 

where  £  denotes  the  realized  sample.  L  is  regarded  as  a  function  of  the 
variance  and  covariance  parameters  in  V  (given  above  in  equations  (2-17)  through 
(2.20))  and  the  similitude  equation  parameters 

Maximum  likelihood  estimates  are  those  values  of  the  parameters  that 
maximize  L  in  such  a  way  that  L  is  guaranteed  positive  not  only  for  the  observed 
values  of  £  but  for  all  possible  realizations  of  the  sample.  In  our  case  this 
means  that  the  maximization  of  L  is  subject  to  the  constraints  that  a  >  0 
and  r  and  T  be  positive  definite.  In  our  notation  the  circumflex  (•» )  will  be 
used  to  denote  the  maximum  likelihood  estimator. 

Rather  than  maximizing  L  directly  one  usually  works  with  the  log-likelihood 

2 

function  X  *  log  L,  which,  upon  inserting  Ha  for  V,  is  found  to  be 

x  -  -  N  log(2ir)  -  £  log  a2  -  I  log  |  H|  -  _L_  (y-Xjj)  'H_1  (£-Xu) .  (3-2) 

2  2  2  2a2 

Parameters  that  maximize  X  also  maximize  L. 


3-3 
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For  Che  present  we  will  use  _9  Co  denote  Che  unknown  variance  and 

covariance  ratios  in  H.  The  functional  dependence  of  X  may  then  be  summarized 
2 

as  \(£,£,cr  ).  Following  Harville  (Reference  5),  we  will  obtain  the 
maximum  likelihood  estimates  j},  a  ,  and  s  in  three  steps: 

(1)  obtain  the  function  ^(.e)  by  maximizing  X  with  respect  to  £  for 

—  2 

an  arbitrary  fixed  value  of  9  (u(0)  does  not  involve  <j  ), 

(2)  obtain  the  ML  estimates  @  and  $  by  the  constrained 

2  -  2 

maximization  of  X*(j9,o  )  =  \(v(e)  ,9^,0  ), 

(3)  obtain  the  ML  estimate  of  u  from 

jl  *  £(£)  *  (3-3) 

•y 

The  same  approach  will  be  taken  to  obtain  REML  estimators  jT,£,  and  a 

2 

except  for  the  use  of  a  restricted  log-likelihood  function  x#(e_,  a  ). 

As  X  is  continuous  and  differentiable  it  can  be  maximized  with  respect  to 

2 

for  fixed  values  of  £  and  o  by  solving  the  system  of  equations 

0.  Using  well  known  rules*  for  differentiation  with  respect  to 
vectors,  we  obtain  the  function  jKs)  as  a  solution  of  the  "normal  equations" 

(X'H"1X)£  -  X'H-1^.  (3-4) 


5See  footnote  5  on  page  1-4. 

*We  use  3 £ ' a / 3£  »  &  and  3£'Gz/3£  *  (G  +  G')£,  where  a  and  £  are 
arbitrary  vectors  and  G  is  a  square  matrix. 
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Since  X  is  Che  N  x  Cp  matrix  of  regressor  values  it  is  assumed  Co  have  full 
column  rank  and  the  matrix  X'H  *X  is  consequently  nonsingular.  The  normal 
equations  then  have  the  unique  solution 

2  -  (X'H_1X)"1X'H'1x-  (3-5) 

Note  here  that  had  we  included  systematic  gage  bias  terms  in  the  model  the 
matrix  X'H  *X  would  have  been  singular.  In  that  event  no  unique  solution  of 
the  full  set  of  similitude  equation  parameters  would  have  been  possible. 

The  close  relationship  between  least  squares  estimators  of  jj_  and  the 
maximum  likelihood  estimator  (assuming  normality  as  expressed  by  (3-D)  is  seen 
by  observing  that  the  same  function  that  maximizes  L  for  fixed  and  o 
also  minimizes  the  weighted  sum  of  squared  residuals  given  by  the  quadratic  form 
(X-XjjVV  ^-X^).  It  follows  from  this  that  equation  (3-5)  may  be  used  to 
represent  a  number  of  different  statistical  estimators  of  ^  obtained  by  least 
squares  and  other  methods.  Some  of  these  are  listed  in  Table  1  and  depend  upon 
the  nature  of  H.  When  H  is  formed  from  the  ML  or  REML  estimates  of  one 
obtains,  as  stated  above,  the  ML  or  REML  estimates  of  ^respectively.  The 
other  estimators  listed  depend  upon  H  having  known  or  prescribed  forms. 

At  present,  similitude  equation  parameters  are  estimated  by  means  of  the 
OLS  (ordinary  least  squares)  estimator,  which  does  not  involve  the  data  variance 
covariance  matrix  H  at  all.  It  can  be  shown  that  all  of  the  other  estimators 
listed  in  Table  1  are  more  accurate  (in  the  sense  of  smaller  variances  or  mean 
squared errors)  than  the  OLS  estimator  when  correlations  exist  among  the  data  and 


ill 
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THE  EXTENDED  NORMAL  EQUATIONS 


In  practice  it  is  usually  computationally  more  attractive  to  find 

formulations  for  £  which  do  not  explicitly  require  the  inversion  of  the  N  x  N 

variance  covariance  ratio  matrix  H  (or  equivalently  V)  as  does  equation  (3-5) . 

There  are  several  ways  by  which  this  can  be  achieved  that  require  the  inversion 

of  matrices  smaller  in  size.  All  of  these  involve  about  the  same  computational 

effort.  A  particularly  useful  formulation  was  published  by  Henderson  et  al. 

(Reference  18)  and  was  investigated  extensively  by  Harville  (see  References  5, 

14,  and  19).  In  addition  to  giving  estimates  of  the  fixed  effects  u,  this 

*  * 

formulation  also  gives  estimates  of  the  means  of  the  random  effects  a  and  1} 

that  are  conditional  upon  £.  These  can  be  regarded  as  estimates  of  the  unknown 

values,  a_  and  j>,  that  are  actually  realized  by  a  and  b  in  the  sample.  We 

will. find  important  uses  for  these  in  the  applications  below.  Furthermore, 

Harville  (Reference  5)  has  shown  how  ML  and  REML  estimates  of  the  variance 

2 

covariance  components  £  and  a  can  be  obtained  from  these  results  with 
little  additional  effort. 


In  order  that  we  might  make  direct  use  of  these  results,  we  rewrite  the 
model  equation  (2-12)  as 


^Henderson,  c.  R. ,  Kempthorne,  0.,  Searle,  S.  R.,  and  Von  Krasigk,  C.  N. , 
"Estimation  of  Environmental  and  Genetic  Trends  from  Records  Subject  to 
Culling,"  Biometrics ,  Vol.  15,  1959,  p.  192. 

5 See  footnote  5  on  page  1-4, 

14See  footnote  14  on  page  3-2. 

^Harville,  D.  A.,  "Extension  of  the  Gauss-Markov  Theorem  to  Include  the 

Estimation  of  Random  Effects,"  The  Annals  of  Statistics,  Vol.  4,  No.  2,  1976, 
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£  ■  X_y  +  Zv*  +  e*, 


and  the  variance  covariance  ratio  matrix  H  given  in  (2-22)  as 


H  ■  I  +  ZDZ ' . 


Here,  quite  obviously,  we  have  defined 


v*  5  (a*',b*')', 


2  5  [U,W], 


■n 


The  work  of  Henderson  et  al.  (from  Searle,  Reference  9,  p.  461)  is  based 

* 

upon  the  joint  likelihood  of  £  and  v  ,  whi-ch  can  be  expressed  as 


g(£,v)  -  gi(z|v)g2(v) 


c q  exp[-  _i_  (i"X£-Zv )  '  (^-Xji-Zv ) ]  expf-  1  v'D“lv], 
2a2  2  a2 


where  cQ  is  a  constant  (function  of  £  and  a  )•  Note  that  although  £  is 

indiscriminant ly  used  in  this  report  to  denote  either  the  random  variable  or 

sample  value  of  the  observation  vector  depending  on  the  context,  our  notation 

with  regard  to  v*  and  v  is  more  explicit.  In  (3-11)  v  indicates  the 

* 

unobservable  realized  value  of  v  .  As  it  is  unknown,  it  may  be  regarded  as  a 

parameter  of  the  model  in  the  same  sense  as  For  fixed  values  of  8  and 
2 

a  ,  maximization  of  (3-11)  with  respect  to  u.  and  v  leads  to  the  system  of 
equations 


9  See  footnote  9  on  page  3-1. 
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where 


x'z 

D“Uz'z 


(3-12) 


(3-13) 


Equations  (3-12)  are  referred  to  as  the  "extended  normal  equations"  or  also 

as  the  "mixed  model  equations."  The  solution  to  (3-12)  gives  the  maximum 

likelihood  estimates  of  £  and  v  for  fixed,  arbitrary  values  of  e  and 
2 

a  .  It  can  be  shown  (e.g.,  see  Serle,  Reference  9,  pp.  459-462)  that  the 

estimator  of  £  from  (3-12)  is  identical  to  that  given  by  (3-5)  and  that  the 

estimator  of  v  is  equivalently  the  ML  estimator  of  E(v*|£).  Hence,  we  can 

denote  solutions  of  (3-12)  as  7  and  7,  and  as  £  and  v  if  the  ML  estimates  of 
•  2 

0  and  o  have  been  used. 


Solution  of  (3-12)  requires  the  inversion  of  the  nonsingular  coefficient 

matrix  A  of  size  Cp+ZMfc+pZJi ,  which  equals  the  number  of  fixed  and 

k  i 

random  levels  in  the  model.  This  will  be  substantially  smaller  than  N.  A 
convenient  way  of  determining  A  ^  is  by  successively  inverting  partitioned 
matrices.  We  define 


(3-14) 


^See  footnote  9  on  page  3-1. 
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Then  by  substitution  of  (3-13)  and  (3-14)  into  AA~ 
Searle,  Reference  20,  p.  210  or  Westlake,  Reference 


*  I  ■  A  *A  (see,  e.g. , 
21,  p.  26)  one  finds 


(3-15) 


A11(x'zB_1) 


(3-16) 


*  Q 


where 


+  z'z 


and 


'  -1 

-  (X  ZB  )  A 


12 


(3-17) 


(3-18) 


(3-19) 


Since  B  and  Q  will  be  used  extensively  in  the  sections  below  they  have  been 
given  special  notations. 

This  method  of  inversion,  thus,  requires  the  inversion  of  B,  the  inversion 

of  D,  and  the  inversion  indicated  in  (3-15).  The  inversion  of  B  will  be 

discussed  below.  The  inversion  of  0  is  easily  obtained  since  it  is  composed  of 

a  diagonal  matrix  and  a  block  diagonal  matrix  consisting  of  repeated  p  x  p 

blocks.  In  practice,  p  will  have  a  value  of  2  or  3  so  that  the  inversion  of  D 

'  *  —i  • 

is  trivial.  The  matrix  X  -  X  ZB  Z  X  in  (3-15)  is  symmetric  and  relatively 
small  in  size  (Cp  x  Cp).  Hence,  it  should  be  possible  to  obtain  A^  rather 

20searle,  S.  R. ,  Matrix  Algebra  for  the  Biological  Sciences  (New  York:  John 
Wiley  &  Sons,  Inc.,  1966). 

21we8tlake,  J.  R. ,  A  Handbook  of  Numerical  Matrix  Inversion  and  Solution  of 
Linear  Equations  (Huntington,  NY:  Robert  E.  Krieger  Publishing  Co.,  1975). 
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rapidly  using  an  approach  such  as  Che  symmetric  Cholesky  method  (see  Westlake, 
Reference  21,  p.  13).  Furthermore,  the  storage  of  and  *  Q  is  aided 
by  the  symmetry  of  both  matrices. 
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b..  -  [u,u+r'1-u'w(w'w+T~1)“1w'u]“1 


(3-24b) 


®12  ’  -B11u'w^w'w+T~1)"1 


(3-23b) 


B„„  -  (I-B*  u'wHw'w+t"1)'1. 


(3-22b) 


One  sees  that  the  only  difficulties  involve  equations  (3-22a)  and  (3-24b); 

-1  1  -1  -1 

the  other  inverse  matrices  required  are  either  diagonal  (r  and  (U  U+r  )  ) 

-1  1  -1  -1 

or  block  diagonal  with  p  x  p  blocks  (T  and  (W  W+T  )  ).  Note  also 

that  r  ^  and  T  ^  are  readily  obtained  from  (or  may  be  used  to  obtain) 

D  The  inversion  indicated  in  (3-22a)  or  (3-24b)  will  be  the  most  time 
consuming  step  in  the  calculation  of  A  Hence,  it  will  probably  be 

important  to.  choose  the  faster  of  the  two  methods  indicated  above.  Since  both 
(3-22a)  and  (3~24b)  are  symmetric,  the  symmetric  Cholesky  approach  is  again  a 
good  choice  of  method. 


3-3  MAXIMUM  LIKELIHOOD  ESTIMATORS  OF  THE  VARIANCE  COMPONENTS 


We  obtain  the  ML  estimates  of  9  and  a  by  maximizing  the  function 
2 

A*(0_,o  ),  which  from  (3-2)  is  written 

\*  •  -  N  log  (2 tt )  -  N  log  a2  -  I  loglHI  -  _L  (y-XjrtV*  -Xu).  (3-25) 

2  2  2  2a2  ~ 

Here  jT  is  given  by  (3-5)  or  obtained  from  the  solution  of  (3-12).  The 
maximization  of  (3-25),  however,  must  be  carried  out  subject  to  constraints  that 
ensure  the  positive  definite  property  of  the  estimated  variance  covariance 
matrices  ?,  T,  and  .  Discussion  of  these  constraints  is  given  below  in 


Section  4-2. 


2 


Using  established  rules*  for  the  differentiation  of  matrices,  we  can  list 

*  .  2 

the  derivatives  of  X  with  respect  to  a  and  an  arbitrary  component 
of  £  as 

1!_  =  _i_  (y-XF)  'H_1(y-XF)  -  -2-  (3-26) 

3<r2  2*4  *  ~  ~  2a2 

il_  -  _L_  (y-XuVH-1  11L  H_1  (y-Xu)  -  I  tr(H_1  2E_).  (3-27) 

3et  2o2  "t  ~  2  36t 

Here,  (3-27)  can  be  most  easily  obtained  by  noting  that 

3X*/39t  *  (3X/3jj)'(3jj/36t)  +  3X/30t  =  3X/3et 

which  is  evaluated  at  u,  *  jl.  This  follows  from  the  requirement  that 
3X/3£  *  0. 

The  maximum  likelihood  estimates  and  *  satisfy,  subject  to  the 

*  2  * 

constraints,  the  "likelihood  equations"  3X  /3o  *  0  and  3X  /36t  *  0 

(for  all  t's),  where  the  derivatives  are  given  by  (3-26)  and  (3-27).  Generally, 
these  equations  are  nonlinear  and  must  be  solved  iteratively.  Recent  numerical 
schemes  for  computing  ML  estimates  of  the  parameters  in  general  mixed  ANOVA 


*We  make  use  of  3log|G|/3s  *  tr(G”^3G/3s)  and  3G"l/3s  *  -G~l  3G/  3s  G-*-, 
where  the  matrix  G  depends  on  the  scalar  s  and  its  trace  tr(G)  is  the  sum  of 
its  diagonal  elements.  Proofs  can  be  obtained  from  Nering  (Reference  22). 

22Nering,  E.  D.,  Linear  Algebra  and  Matrix  Theory  (New  York:  John  Wiley  & 
Sons,  2nd  Ed.  1970). 
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models  have  employed  the  Newton-Raphson  (N-R)  procedure  (see  References  11  and 

15)  or  a  combination  of  the  N-R  procedure  and  Fisher's  method  of  scoring 

(Jenrich  and  Sampson,  Reference  13).  The  method  of  scoring  is  a  modification  of 

the  N-R  algorithm  which  uses  the  matrix  of  expected  values  of  the  second 

derivatives  in  place  of  the  matrix  of  second  derivatives  (the  Hessian).  The  N-R 

method  is  the  more  efficient  approach  and  the  method  of  choice  when  the 
.  .  .  ★ 

log-likelihood  function  X  is  approximately  quadratic,  but  Jennrich  and 
Sampson  suggest  that  it  be  backed  up  by  the  method  of  scoring  under  poor 
starting  conditions  or  if  the  N-R  iterates  begin  to  diverge.  We  discuss  these 
techniques  more  fully  in  Chapter  4. 

Employment  of  these  iterative  schemes,  thus,  requires  the  availability  of 
second  derivatives  of  X*  and  the  expected  values  of  second  derivatives. 

Methods  for  computing  these  quantities  have  been  developed  by  several 
researchers  (References  10,  11,  12,  and  14),  but  the  most  useful  results  have 
been  obtained  by  Harville  (Reference  5)  who  has  shown  how  they  may  be  extracted 
from  the  results  of  the  previous  section.  Hence,  the  development  below  is  based 
largely  upon  Harvi lie's  work. 


11-See 

footnote 

11 

on 

page 

3-1. 

l^See 

footnote 

15 

on 

page 

3-2. 

13  See 

footnote 

13 

on 

page 

3-2. 

l°See 

footnote 

10 

on 

page 

3-1. 

l^See 

footnote 

12 

on 

page 

3-1. 

l^See 

footnote 

14 

on 

page 

3-2. 

^See 

footnote 

5  or 

i  page 

1-4. 
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•ft 

The  second  derivatives  of  X  can  be  obtained  from  (3-26)  and  (3-27)  and 
listed  as  follows: 


3  x  -  =  -  _JL  (Z-Xu)  ’H-ki-Xu)  +  JL 

3a2  3o2  c6  2 


(3-28) 


3eg3a2 


_L_  (v-XID'h-1  [is_  (l-2PH)l  H-1  (y-Xu) 

2ct^  [3«s  J 


3  9g  3  9 


_  =  _L_  (y-XiT)  'h-1  [&-_  -  2  i2_  P  iS_l  H-1  (y-Xu) 
t  2a2  L39s 39t  39s  30tJ 


where 


-1  -l  '  -1  -1  1  -1 

P  5  H  -  H  X(X  H  X)  X  H  . 


(3-29) 


!  -  I  tr 

{h-1  f  32H  -  3H  H-l 

3H  I/ 

9 

(3-30) 

2 

{  [3es39t  39s 

30tJ( 

(3-31) 


In  deriving  these  expressions  use  is  made  of  the  result  3P/39  * 

s 

-P(3H/30g)P,  which  can  be  shown  by  differentiating  (3-31)  and  performing 
some  simple  algebra.  It  should  also  be  noted  that  P£  *  H  *(£~XjT)  and  that  P 
is  symmetric. 


To  obtain  the  expected  values  of  these  derivatives  one  requires  the 
expression  for  the  mean  of  a  quadratic  form  (see  Searle,  Reference  9,  p.  55) 

E^'gj.)  *  tr(Gv)  +  uVGXu.,  (3-32) 

which  is  true  for  a  general  matrix  G.  Employing  (3-32)  and  the  facts  that 


^See  footnote  9  on  page  3-1. 
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(PH;  -  PH, 


(3-33) 


PX  *  0, 


(3-34) 


and 


tr(FG)  =  tr(GF) 


(3-34) 


for  arbitrary  matrices  F  and  G,  one  finds  that 


s  ( J2i _ ^  -  -  JL  (N-cp)  +  _JL 

\3<r23a2/  cA  2ct4 

/A*  V  -  _J_  tr  (?  iS_\ 
\39g3s2/  2a2  \  39s/ 

V  I  tr  jp  r 

\3es3et/  2  (  L 


-  I  tr 
2 


3  H  _  2  3JL  P  1S_ 
90g86^  30g  30^ 


fH-l  f  32H  _  3H_  H-1  3H_ 

I  L30s3®t  30s  30t. 


(3-36) 


(3-37) 


(3-38) 


The  derivatives  in  equations  (3-26)  through  (3-30)  and  the  expectations 
given  by  (3-36),  (3-37),  and  (3-38)  can  be  simplified  by  using  certain  results 
that  we  now  indicate  (see  Harville,  Reference  5).  Many  of  these  results  derive 
from  the  use  of  the  matrix  identity 


(E+FG)"1  =  e”1-E-1F(I+GE"1F)"1GE~1 


(3-39) 


^See  footnote  5  on  page  1-4. 
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(obtained  by  representing  (I+FGE  *)  ^  by  a  geometric  series,  see  Bartlett, 
Reference  23).  Using  (3-7)  we  find 


1  -1  •  _i  • 

ZH  *  (I+Z  ZD)  Z 


(3-40) 


»  •  _i  • 

Z  P  *  (I+Z  SZD)  Z  S, 


(3-41) 


where  S  =  I-X(X  X)"  X  . 


(3-42) 


From  (3-42)  it  follows  that  SX  =  0.  Next  we  find  several  expressions  for 


-  „-l- 

v  =  D  v  , 


(3-43) 


which  is  a  quantity  introduced  by  Harville  to  facilitate  the  development  of  the 
\*  derivatives.  Using  the  above  relations  with  the  result  from  (3-12)  and 

_f  t  _  « 

(3-14)  that  v  *  (A,.X  +  A_„Z  )£  one  can  determine  that 


'  -1  ' 

2  ■  (I+Z  ZD)  Z  (£-Xjp 


z  h  ^-xjr) 


(3-44) 


(3-45) 


(i+z'szdTVsx 


(3-46) 


^Bartlett,  M.  S.,  "An  Inverse  Matrix  Adjustment  Arising  in  Discriminant 
Analysis,"  Annals  of  Mathematical  Statistics,  Vol.  22,  1951,  p.  107. 
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Finally,  from  (3-39)  and  (3-44)  we  obtain  the  useful  result 
H  * (Z-Xu)  -  ^-XF-Zv  =  ¥  , 


(3-47) 


which  is  the  estimated  error  vector  or  vector  of  residuals,  and  from  (3-5), 
(3-31),  and  (3-34)  we  find 


(y-XjpV^r-XjD  *  x'pZ  • 


(3-48) 


With  these  results  and  remembering  that  H  depends  on  0  by  way  of  D,  we 
can  list  the  sought  after  derivatives  as 


3  A 


3o2  3o^ 


1  y '  e  -  N 
“  3a2 


(3-49) 


I 


3  A 
30t  2 


_L  v'  12-  v-Itr  /  [  I+Z '  ZD]  _1Z  ’  Z  £jLl 
j„2  ~  30t  ~  2  {  ^0^  ( 


30t  ) 


2  * 

3  A 

3a2  3<j2 


-  _L  y'e  +  JL 
.  06  2o^ 


(3-50) 


(3-51) 


2  * 

3  A  .  _1  v*  _3D_  -y 


39t3o2  2  30t 


(3-52) 


2  * 

3  A 

39s30c 


-  JL  v'  i2_  [I+z'sZDj-iz'sZ  12_  V 


t2  “  39, 


30t 


1  tr  |[i+z'zd]"1  z'z  i£_  [i+z’zd]-1  z'z  i£_\ 

2  <  39s  30tf 


(3-53) 
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(  *Z  X  \  -  -  1_  (N-Cp)  +  JL 
\3  o23o2y  a^  2  a ^ 

(JL±J\  *  -  _1_  Cr  /  [I+z'sZDJ^z'sZ  12-1 
3 0t3o2/  2a2  ^  3®t  f 

*  -tr  ([I+Z'SZDJ^Z'SZ  12-  [i+z'sZDj-iz'sZ  ?/>  I 
3®s 3®t y  ®®s  3®t > 

+  I  tr  /[i+z'zd]"1z'z  12_  [i+z'zd]  z'z  12_l  . 

2  (  36g  30tJ 


(3-54) 


(3-55) 


(3-56) 


It  remains  to  express  these  derivatives  in  terms  of  quantities  derived  in  the 

previous  section.  Those  derivatives  that  appear  complicated  depend  upon 

matrices  of  the  forms.  [I+z'sZDj^z'sz  ^  and  [I+z'zdJ'^z'z  .  The 

38fc  30t 

latter  of  these. is  readily  seen  to  depend  on  the  matrix  B  a  calcula- 


tional  scheme  for  which  was  indicated  in  equations  (3-22)  through  (3-24).  Hence, 
from  (3-18)  we  can  write 

(i+z'ZDj-J-z'z  12-  *  D^B^z'z  12-  .  (3-57) 

30t  30t 

The  former  of  these  matrices  has  been  shown  by  Harville  (Reference  5,  p.  326)  to 
depend  upon  the  matrix  Q  given  by  equation  (3-19).  It  can  be  shown  by  expansion 
of  (3-19)  and  the  use  of  (3-39)  and  (3-42)  that 

[1+z'sZD]-1  =  D_1Q,  (3-58) 

^See  footnote  5  on  page  1-4. 
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and,  therefore,  that 

[i+z’szd]"1z'sz  UL  »  d-1oz'sz  i£_  . 

30t  39t 

Now,  by  means  of  the  identities 


(3-59) 


b'z'z  5  -f B~1D~1-I ]  (3-60) 

and  Qz'sz  =  -[QD"1-!]  ,  (3-61) 

obtained  from  I  ;  [D  *+Z  SZ]  *  and  I  =  [D  *+Z  Z]  *B,  it  is 
possible  to  put  (3-57)  and  (3-59)  into  very  similar  forms.  We  find 

[I+Z'ZDJ-J-Z'Z  i£_  «  -[D-J-B-i-ljD-1  iL. 

30t  30t 

and  II+Z'SZDI-VSZ  jUL.  -  IS-  . 

30t  30t 

As  a  consequence  of  this  similarity  of  form  we  now  need  only  continue  the 
development  of  (3-62),  say,  and  apply  the  results  to  (3-63)  by  substituting  Q 
for  B  ^ . 


(3-62) 

(3-63) 


Equation  (3-62)  can  be  simplified  by  examining  the  structural  details  of 
its  matrices  and  partitioning  them  according  to  both  gage  classes  and  explosive 
classes.  For  D  *  Block  Diag  (r,T)  we  define 


T.  s  Diag 


(rk»*-*.Yk) 


Vm. 


(3-64) 


3-20 
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and 


T.  =  Block  Diag 


(3-65) 


Hence,  from  (2-17)  and  (2-19)  the  stated  partitioning  of  D  may  be  expressed  by 


r  =  Block  Diag  ( , Tj, • • • , rR) 


(3-66) 


and 


T  i  Block  Diag  (T^ ,T2> . . . ,TC)  . 


(3-67) 


In  a  like  manner  we  now  partition  B-*'.  From  (3-20)  and  3-21)  we  have 


u'u+r*1 

u'w  ■ 

-1 

'Bll 

b12  " 

.  w'u 

W'W+T-1 

mi 

-B21 

B22_ 

(3-68) 


where  it  is  convenient  to  write  B^  for  B^ •  Then  the  above  partitioning 

•  . 

may  be  indicated  by  writing 


B11  * 

^ll(l.l) 

^1 1  (  2 , 1 ) 

*n(l,2) 

^1 1 ( 2 , 2 ) 

•  •  • 

•  •  • 

BU(1,K) 

1U(2,K) 

9 

(3-69) 

. 

BU(K,1) 

Bn(K,2) 

•  •  • 

Ill(K,K) 

B22  * 

B22(l,l) 
1^22(2  >  i  ) 

B22(l,2) 
^22^  >  2) 

•  •  • 

•  •  • 

B22(1,C) 
TT22(2  ,C) 

9 

(3-70) 

B22(C,1) 

122^,2) 

•  •  • 

B22(C,C) 

b12  “ 

ri12a.i) 

B12(2,l) 

Bi2(1,2) 
5*12  ( 2 ,2) 

•  •  • 

•  •  • 

Ul2(l ,0 

7j2(2,C) ! 

9 

(3-71) 

_Bi2(K,1) 

Bi2(K,2) 

•  •  • 

Bi2(K,C)_ 

a  similar  structure  for 

®21  *  *12’ 

Note  that  from  the  symmetry  of 

we  have  B^^(k^,k2)  ■  (B 

11(k2,k1)) 

’  B22 

_  I 

(i1,i2)  =  (B22(i2’il))  *  and 

NSWC  TR  82-74 


I 

I 


N 


4 


At  this  point  it  is  useful  to  distinguish  between  the  different  parameters 

in  _0.  These  were  indicated  previously  in  equations  (2-18)  and  (2-20)  as  the 

gage  class  calibration  error  variance  ratios  y^,k  *  1,...,K  and  the 

explosive  class  performance  variance  covariance  ratios  in©,,i  =  1,...,C.  Let 

*  1 

us  now  indicate  the  tth  unknown  element  of  ©^  (and  also  T^)  as  t^. 

Here,  the  actual  assignment  plan  is  arbitrary.  There  will  be  p!  of  these 
parameters  for  each  explosive  class. 


Choosing  0t  *  r- t ,  the  derivative  3D/36t  =  3D/3t£t  in 
(3-62)  is  given  by 


3D _ 

3rit 


0 

» 

3T/3rit_ 


(3-72) 


where 


3T . 

aT  -  *  Block  Diag  (0,...,0,  .  *■  ,Q,...,0)  . 

3rit  3tit 


Consequently, 


(3-73) 


D-l  jD  -  Block  Diag  (0,...,0,T71  _ L  ,0,...,0)  .  (3-74) 

®Tit  *•  3Tit 

Also,  the  matrix  [D  *B  *  -  I]  in  (3-62)  consists  of  the  rows  of  the 
partitioned  B  *  matrix  premultiplied  in  succession  by  the  diagonal  blocks  of 
D~1  with  identity  matrices  of  the  forms  and  Ij^  subtracted  along  the 


diagonal.  Postmultiplication  of  this  by  (3-74),  then,  gives  a  matrix  consisting 
of  a  single  column  of  submatrices  as  shown  in  (3-75).  And  by  a  similar  argument 
one  obtains  (3-76).  In  these  the  dashed  lines  indicate  the  partitions  between 


*Tit 
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explosive  classes  and  gage  classes.  As  stated  above,  these  equations  apply  to 
equation  (3-63)  by  replacing  B  ^  or  B  matrices  by  similar  Q  terms. 


These  results  can  now  be  used  to  simplify  the  derivatives  in  (3-49)  through 
(3-56).  In  doing  so  we  use  the  definition 


(3-77 


where  £  corresponds  to  the  gage  class  elements  and  n  to  the  explosive  class 
elements.  Furthermore,  we  choose  T  and  to  be  partitioned  in  the  manner 

_  _ i  i  _ i  t  it  i  » 

discussed  above  so  that  £  =  (^  ,£2 , . . .  .j^)  and  ][  =  (jn,i  i]!? » •  •  •  •Sp) »  where 
^  and  correspond  to  the  kth  gage  class  and  the  ith  explosive  class 
respectively.  The  derivatives  then  become 


*  y'e  -  _JL 
2ct4  2a2 


*  '  /  ^  M 

3*  *  1  EE  +  tr  <  if  (k,k) >  -  _ii 

3Yk  2  "2  **  ■»  .2  'll  ’  y 


ii ■  iL-  r\'  —L.  rf  +  I  tr  ITt-1  B  (i,i)  -  Ij  1  T_1  -!Il 

3tit  2o2  1  3xit  i  2  (L  i  22  ij  i  3rit 


(3-78 


(3-79 


(3-80 
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J3T . 

(T_1  0  Ci , i >  “  Ij  )  T*1  _ L 

i  22  i  i  3ris 


3T .  ) 

(T-1  0  (i,i)  -  It  )  T"1  — 

i  22  i  i  3r£t  ) 


(  3T. 

+  I  tr  < (T-1  B  (i,i)  -  I j  )  T-1  _ L_ 

2  /  i  22  i  i  3t£s 


(T-1  B  (i,i)  -  It  )  T"1  !TL  t 
i  22  i  i  3-rit  ) 


(  3T 

-  tr  )t-1  Q  (i,h)  T"1  h 
/  i  22  h  3t^s 


T~1 

Q  (h ,  i )  T'1 

3Ti  j 

h 

22  i 

3*it! 

+  I 

tr  )t-1  B  ( 

i/h)  T*1 

3T. 

.  n  • 

2 

(  i  22* 

h 

3rhs 

T-1 

B  (h,i)  T”1 

!!i_j  , 

h  *  i 

h 

22  i 

3Ut) 

(3-95 


(3-96 


The  final  simplification  of  the  derivatives  is  achieved  by  exploiting  the 

repetition  of  ©^  within  T^.  As  indicated  in  (2-20)  and  (3-65)  ©^  is  a 

p  x  p  matrix  which  is  entered  along  the  diagonal  of  for  each  of  the 

shots  of  the  explosive  class.  To  take  advantage  of  this  structure,  then,  we 

further  partition  V  and  the  matrices  B  ^  and  0  according  to  individual 

shots.  We  will  denote  this  by  subscripting  the  explosive  class  index  or  indices 

of  the  particular  submatrices.  For  example,  Bj^k.i^)  will  refer  to  the  p 

columns  of  B10(k,i)  that  pertain  to  the  uth  shot,  and  "B„„(h  ,i  )  will 
12  ?2  u  v 

denote  the  p  x  p  submatrix  of  B^jfh.i)  that  pertains  to  the  uth  shot  of  the 
hth  explosive  class  and  the  vth  shot  of  the  ith  explosive  class. 
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Employing  this  notation  we  can  express  the  derivatives  in  forms  highly 


suited 

for  ' 

* 

3X  = 

1  I 

3  a2 

2o* 

★ 

3  A  = 

1  I 

3Yk 

2a2  ^ 

* 

3A 

1 

3xit 

2a2 

+  I 

2 

0  * 

3  X 

3  a2  3a2 

2  * 

3  A 

3Yk3<*2 

2  * 

3  A 

'e  -  JL 

2<j2 


,  30. 


+  I  Y!  tr  I  eT1  ®  (i  ,i  )  -  I  1  oT1  — i-i 

2U  (l  l  22  uu  pji  I 


4z« 


2  tA 


-  -  I  2  5!  y—  I. 

3xit3a2  ^  u  lu  3xit  Lu 


(3-97) 


(3-98) 


(3-99) 


(3-100) 


(3-101) 


(3-102) 


3Yk3Yk 


_JL_  t  fJL  0  <k,k)  -  lu  ]  K 

2^  Tc  l_Yk  n  \l  "k 


+  _i_  tr  )r_l  B  (k,k)  -  Im  "ft 


(3-103) 
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3  A  -  =  -  i _  s'  Q  (k,i)  _£ 

37k3Yi  2ct2y  Y  *  H  * 

k  l 


+  1  tr  /b  (£,k)  B  (k,jO>  ,  k  *  l 

11  ' 


a  x  _  =  1-  £  t  Q  (k,i  )  ©71  — i-  w 

3Yk3xit  2a2Yk  u  1^  u  i  3t£t  iu 


3tj  a  3r 


is  aTit 


„1_  V  tr  •  ©“1  B  (i  ,k)  B  (k,i  )  ©~1  — L_ 
2yfc  u  i  21  u  12  u  i  3fit 


=  -L-T,  £  n'  Te-l  Q  (i  ,i)-lle-l^-n 

2a2  u  v  _iu  3xis  Li  22  u  v  pj  l  3rit  iv 

+  i  2  tr  IT©?1  ®  (i  ,i  )  -  I  1  ©71  — i- 
2  u  (I  1  22  uu  pji  Sxis 

fe“l  B  (i  ,i  )  -  I  1  ©7I  — i—  \ 

L  1  22  u  u  PJ  1  3tit  \ 

=  J_  T  T  n  _JL_  ©'I  Q  (h  ,  i  )  e-1  _1_  n. 

o-2  77  t-  ~h„  3tv, q  h  22  u  v  1  3rit  “Yv 


3xhs3Tit  2 u  v  ~~ 3xhs 


+  i  V  tr  <0-l  B  (i  ,h  )  ©“I 
2  u  (i  22  uu  h  3ths 


0*1  B  (h  ,  i  )  0-1  — i—  >  ,  h  *  i 

i  22  u  u  i  3T£t  ) 


(  32*  \  *  -  L_  (N-Cp)  +  JL 

\3a23a2y  2a^ 

/  3  X  \  =  1  Q  tr  |o  (k,k)|  -M^l 

\3Yk3c2,/  2a2Yk  |rk  <  “  '  J 


(  -104) 


(3-105) 


(3-106) 


(3-107) 


(3-108) 


(3-109) 
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(2  * 

■  UL 

3-ti  f  3c 


-]  a  — L_  £  tr  jf©?1  Q  (i  ,i  )  -  I  I  eT1  — i-l 
7  2a2  u  jL  1  22  u  u  pj  i  3tit  ( 


f  3  V  ]  =  -  i_  tr  Q  (k 

yrkt-rkj  ri  )j_7k  11 


»k)  -  Iv 


_  2  \ 

+  _L  tr  JQ  B  (k,k)  -  IM  1 
2y£  |Lrk  11  M*J  \ 


/  3  -  *  \  *  -  _i _  tr  i  Q  ( t 

y3Yk3Yjl  J  yjjyJ  ^  11 


,k)  Q^Ck.^j 


2y£y$ 


tr  s B  ( i,k)  B  (k 


,i)|  , 


k  *  l 


(  t2*  )  =  -  i_  V  tr  le?!  Q  (i  ,k)  Q  (k,i  )  ©7! 

y3Yk3Tit J  Y^  u  (  1  21  u  12  u  i  3rit 

+  _1_  tr  I®?1  ®  »k)  B  (k,i  )  ©7I  _ 

2  u  (  1  21  u  12  u  i  c 

(  J\  =  -Itr  jf©:l  Q  (i  ,i  )  -  I  1  ©71  !!L 

\3Tis 3Tit/  u  |[  1  22  u  u  PJ  1  3ris 

[©7I  Q  (i  ,i  )  -  I  1  ©7I  ffi-l 

I  1  22  u  u  pj  1  3T£t  | 

*  1  Z  tr  )[©-!  F  (i  ,i  )  -  I  1  ©-1  !!l 
2  u  (L  i  22  u  u  pj  1  3fi£ 

|©-1  B  (i  ,i  )  -  I  1  ©_1  — i — ( 

(_  1  22  u  u  pj  x  3xit  , 


1  30i 


(3-110) 


(3-111) 


(3-112) 


(3-113) 


(3-114) 
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(-&L) 

\3Ths  3Tit/ 


£  tr  I©”1  Q  (i  .h  ) 

U  i  22  u  u 


1  1  i  3@. 

eh  ar-eh  Q22(hu*V  e'i  rr 

hs  it 

+  1  j®il522(Vh«> 


,3©,  ,  _  ,  a©.  ) 

eh  ar-  *1  B22(hu*V  •;  sr;  •  h  * 1 

hs  it  ) 


(3-115) 


In  these  expressions  the  sums  on  u  and  v  range  from  1  to  the  number  of  shots  in 
the  explosive  class  (J^  for  the  ith  class). 


3-4  RESTRICTED  MAXIMUM  LIKELIHOOD  ESTIMATORS  OF  THE  VARIANCE  COMPONENTS 


It  is  well  known  that  maximum  likelihood  estimators  of  variance  components, 

symbolized  above  by  a  and  J^,  are  biased,  i.e.,  E (£z )  *  a  and  E(j})  *  e. 

Hence,  the  distributions  of  the  estimators  are  not  centered,  in  the  sense  of  the 
means,  about  the  true  values  of  the  quantities  being  estimated.  For  small  data 
samples  the  bias  can  lead  to  substantial  errors  in  estimation.  Corbeil  and 
Searle  (Reference  15),  for  example,  show  for  a  sample  of  size  16  that  ML 
estimators  can  underpredict  estimators  that  are  known  to  be  unbiased  by  a  factor 
of  two.  The  restricted  maximum  likelihood  (REML)  method  is  an  attempt  to 
overcome  this  problem  of  bias.  It  has  been  demonstrated  to  produce  unbiased 
estimates  when  applied  to  a  number  of  different  linear  models  (see  Harville, 
Reference  5),  and  the  property  of  unbiasedness  is  believed  by  some  to  be  a 


^See  footnote  15  on  page  3-2, 
^See  footnote  5  on  page  1-4. 
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general  property  of  the  method.  Corbeil  and  Searle  (Reference  16),  however, 
have  shown  that  this  may  come  at  the  expense  of  estimator  efficiency.  Because 
of  the  unbiasedness  property  of  the  REML  variance  components  estimators  and 
their  close  association  with  the  ML  estimators,  we  include  them  in  this  report. 

2 

The  REML  estimators  of  o  and  e_  are  based  upon  the  likelihood 
function  associated  with  certain  linearly  independent  combinations  of  the  data 
which  possess  zero  expectations.  Such  sums  of  the  observations,  of  which  there 
is  a  total  of  N-Cp  in  the  present  case,  are  known  as  the  error  contrasts.  It  is 
argued  (Patterson  and  Thompson,  Reference  17)  that  they  may  be  thought  of  as 
containing  all  of  the  variance  component  information  and  should,  therefore,  form 
the  basis  of  estimation. 

If  we  let  R£  denote  a  particular  set  of  the  error  contrasts,  it  has  been 
indicated  by  Harville  (Reference  5)  that  to  within  an  additive  constant  the 
log-likelihood  function  for  R£  may  be  written  as 

X*  .*  -  (N-Cp)  iog  a2  -  I  log|  H{  -  I  log)  X  ,H”1Xj  -  _L.  (£-XjT)  V1  (£-X][) .  (3-116) 


l^See  footnote  16  on  page  3-2. 

*2See  footnote  17  on  page  3-2. 

*Harville  cites  his  1974  paper  (Reference  24)  for  this  derivation,  but  it  is 
derived  in  a  Bayesian  context.  However,  the  author  has  been  able  to  verify 
this  expression  in  a  classical  setting  using  rather  straightforward  variable 
transformation  theory  along  with  several  matrix  relations  published  by  Harville 
in  his  1974  paper. 

2^Harville,  D.  A.,  "Baye  sian  Inference  for  Variance  Components  Using  Only 
Error  Contrasts,"  Biometrika,  Vol.  61,  1974,  p.  38. 

^See  footnote  5  on  page  1-4. 


2 

The  REML  estimates  of  9.  and  a  are  those  values  that  maximize  (3-116).  We 
now  show  how  the  derivatives  required  for  the  maximization  of  X  by  the  N-R 
and  Fisher  scoring  optimization  schemes  may  be  obtained  from  the  derivatives  of 


In  the  same  manner  that  we  obtained  the  X  derivatives  given  by  (3-26) 
through  (3-30),  and  (3-36)  through  (3-38)  we  can  write 


3  a2 

2  CT4 

3X#  , 

1 

38^ 

2  or2 

»V 

3a2  3a2 


- *  -  i_  (£-Xl[)  V^-XjT)  +  -1-  (N-Cp) 


(3-117) 


(3-118)* 


(3-119) 


30s3o2 


_L_  (y-Xu")  'h-1  fiS-  (I-2PH)]  H_1  (y-Xy ) 

2a4  I39*  J 


30s30t 


«  _J_  (y-XjT)  'h-1  L  32H  -  2  ijL  P  JjL  H_1(y-Xu) 

«t  2  a2  13  9  3  3  0 1  30s  39t 


(3-120) 


-I  trip  r»2H  -  3H_  P  3_H  1) 

2  (  I 39s30t  30s  30tJj 

E  ■  -  i_  (N-Cp)  +  (N-CP_)  *  -  JL  (N-Cp) 

y3o23 a2/  a4  2a4  2a4 

♦Note  that  if  x#  5  g*(o2,  9.,  £(0.))  and  g  -  g#(o2,  _9,  u.) ,  that  £  is  a 
solution  to  3g/3u  *  0.  Hence,  3X^/30f  *  3g/3©t  evaluated  at  a  *  £. 


(3-121) 


(3-122) 
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l30»3o2 


»  -  - 1  tr  |p  iiL 
/  2a2  /  3»3 


(3-123) 


(jV.1  .  1  tr  |  p  f  32H  .  2  ILp  3H  ~R 

3es3  0ty  2  (  l_3e830t  30g  30tJj 

-  I  tr  (?  32H  -  3H_  H“1  3H_1> 

2  (  L30S39t  3  ®3  38tJ) 


—  i  tr  |  P  3H  p  3H 
2  I  30s  30t 


(3-124) 


It  is  observed  that  equations  (3-117)  through  (3-124)  differ  only  slightly 
from  equations  (3-26)  through  (3-30)  and  (3-36)  through  (3-38).  The 
expectations  of  derivatives  are  shown  above  both  in  a  way  that  emphasizes  this 
similarity  and  also  in  their  algebraically  simplified  forms.  Close  inspection 
of  these  results  reveals  that  we  can  immediately  write  down  the  REML  equivalents 
to  equations  (3-97)  through  (3-115)  by  simply  using  the  submatrices  of  Q  where 
those  of  B  1  appear  and  by  employing  (N-Cp)  in  certain  places  where  only  N 
appears.  By  proceeding  in  this  manner  and  simplifying  the  results  we  obtain  the 
computational  forms  of  the  derivatives  as 


li_  -  y'e  -  JL  (N-Cp) 
3o2  2o*  2a2 


(3-125) 


*  1  c'T  ♦  1  tri<)  (k,k)j  - 

3Yk  2"2  Kk  jvZ  (  11  )  2yv 


(3-126) 


llL  -  _L  £  I.'  n 

®Tit  2a2  u  *u  3Tit  *u 

+  I  23  tr  fcrl  0  (i  ,i  )  -  I  1  e_1  . \ 

2  u  (L  1  22  u  u  pj  i  3rit  \ 


(3-127) 


l 
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3o23o2 


2  # 

3  X 

3Tit 3cr2 
.2.# 


I_  y +  _1_  (N-Cp) 
c*  2  a1* 


i  ft 


1  ~ 

1  L  I  — = —  n. 

2  u  iu  3t£t  1u 
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CHAPTER  4 
NUMERICAL  METHODS 

Following  Jennrich  and  Sampson  (Reference  13)  we  propose  the  use  of  the 
combined  Newton-Raphson  and  Fisher  scoring  procedures  to  obtain  the  ML  or  REML 
estimates  of  the  variance  compo*ents.  Satisfaction  of  the  constraints  on  the 
iterated  variance  covariance  matrices  will  be  handled  by  the  interior  penalty' 
function  technique.  While  the  expressions  below  may  be  used  interchangeably  to 
obtain  either  ML  or  REML  estimates,  they  will,  for  the  sake  of  brevity,  be  given 
in  most  cases  in  terms  of  the  ML  notation  only. 

4-1  UNCONSTRAINED  OPTIMIZATION 


.  Concise  descriptions  of  the  Newton-Raphson  and  Fisher  scoring  algorithms 
may  be  stated  in  terms  of  the  parameter  vector 

m.  =  ( » Y1 » •  •  • » Yk* T11  ♦ « •  •  >  T , » •  • »  Tqj  , » .  • ,  ,  (4—1) 

/ 

where  ir  5  pi.  The  N-R  method  is  an  iterative  procedure  that  corrects  an 
initial  guess  or  the  value  of  the  iterate  after  i  steps,  ,  by  an 
amount 


4(i)i 


I 


^See  footnote  13  on  page  3-2. 


(4-2) 


4-1 
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where  Am^  *  £Li  +  l  21?  t*ie  gradient  vector 


(*  *  *  * 

3X  3X  ...  3X  ax 

2a2  3n  3yk  3ni 


*  *  y 

,ii_ . 1L_] 

3tC1  3  TC  itj 


(4-3) 


and  3C  is  the  Hessian  matrix  defined  by  equation  (4-4)  shown  on  the  opposite  page. 
Note  that  in  (4-2)  VX*  and  JCare  evaluated  at  w. .  The  method  of  scoring  is 
identical  to  the  N-R  algorithm  except  that  EQf©  is  used  in  place  of  JC. 


Both  the  N-R  procedure  and  the  method  of  scoring  approximate  the  X 
function  at  each  step  of  the  iteration  by  a  quadratic  function.  In  the  case  of 
the  N-R  method  this  is 


X*(w  )  “  X*(w  )  +  (VX*)  Aw  +  i  (Aw  ) '  Aw  , 
“i+1  i  l  i  2  i  i  i 


(4-5) 


which  is  recognized  as  the  Taylor  expansion  of  X*  about  the  point  up  to 
the  second  order  term.  Upon  taking  the  gradient  of  X*(w^+^)  we  get* 


VX*U  )  =  VX*  +3QAm 
i+1  i  l  i 


(4-6) 


which  must  be  zero  at  a  maximum.  Hence,  equating  (4-6)  to  zero  and  solving  for 
Am.  gives  (4-2).  In  order  for  (4-5)  to  approximate  X  near  a  maximum, 
as  opposed  to  a  minimum  or  saddle  point,  the  Hessian  must  be  negative 

I 

definite,  i.e.,  zJGz  <  0  for  any  vector  z  ^  0  (in  particular  Am.),  since 

★ 

any  excursion  away  from  the  maximum  must  result  in  a  decrease  in  X  . 


♦Use  footnote  on  page  3-4  and  symmetry  ofJC* 
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The  fact  that  away  from  the  maximum  JCneed  not  be  negative  definite  can 

cause  the  N-R  interate  to  diverge  from  the  maximum  if  the  starting  conditions 

are  poor.  It  is  for  this  reason  that  Jennrich  and  Sampson  (Reference  13)  employ 

the  method  of  scoring  in  the  initial  step  of  the  iteration  and  whenever  the 

process  appears  to  diverge.  Since  E(JO  is  nonpositive  definite  (<0)  and  in 

our  case  will  most  likely  be  negative  definite  (<0),  the  scoring  step  Au^.  = 

—  1  * 

-[ECJO]  VX^  will  at  least  locally  always  be  in  the  direction  of 
.  .  * 

increasing  X  .  To  show  this  we  examine  the  component  of  A  m  the 
direction  of  VX£ .  Using  Figure  6,  this  is 


|Aw.jcos<j>  *  (Aw.  )  '  VX_*/|  AX_*|  =  -(VX*)  '  [ECJO  ]-1 


(4-7) 


Hence,  cos  <(>  >  0  and  Au>£  has  a  component  in  the  direction  of  increasing  x  . 


Figure  6.  X*  contours  and  parameter  vectors  associated  with  the  Newton-Raphson 
and  method  of  scoring  optimization  schemes 


^See  footnote  13  on  page  3-2. 
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To  prevent  the  step  Atd^  from  overshooting  the  local  region  of 
* 

increasing  X  the  step  is  often  written  in  the  modified  form 

Awi  =  -p£ [E(J©  ] -1  ,  (4-8) 


where  0  <  p-  <  1.  The  value  of  p^  can  be  chosen  in  various  ways  (e.g. 

see  Fox,  Reference  25)  but  Jennrich  and  Sampson  (Reference  13)  simply  halve  the 

★ 

preceeding  step  until  an  increased  value  of  X  is  obtained. 


4-2  CONSTRAINED  OPTIMIZATION 


It  is  expected  that  for  most  sets  of  data  the  unconstrained  N-R  and  Fisher 
scoring  algorithms  will  produce  estimates  of  that  lie  in  the  feasible  region 
of  jo  space  —  that  is,  the  region  corresponding  to  positive  definite  variance 
covariance  matrices.  When  such  is  not  the  case,  however,  we  need  a  method  by 

it 

which  the  maximum  of  X  within  or  at  the  boundary  of  the  feasible  region  can 
be  found.  A  popular  and  general  method  for  effecting  such  a  solution  is 
referred  to  as  the  interior  penalty  function  technique  (see  e.g.  Reference  25 
and  26).  We  will  use  a  version  of  this  method  proposed  by  Carroll  (Reference 

it  it 

27)  that  maximizes  the  function  A  (w)  instead  of  X  ,  where  A  is 
defined  as 


25gox,  R.  L. ,  Optimization  Methods  for  Engineering  Design  (Reading,  MA: 
Addison-Wesley  Publishing  Co.,  1973). 

l^See  footnote  13  on  page  3-2. 

^Aoki,  M. ,  Introduction  to  Optimization  Techniques  (New  York:  The  Macmillan 
Co.,  1971). 

2?Carroll,  C.  W. ,  "The  Created  Response  Surface  Technique  for  Optimizing 
Nonlinear,  Restrained  Systems,"  Operations  Research,  Vol.  9,  169. 
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A*((u)  *  X*(u)_)  -  ?(w) 

and  ?(w)  =  ^  • 

l 

Here  £_  =  {te^}  is  a  set  of  positive  constants  and  we  assume  that  the  set  of 

constraint  relationships  can  be  expressed  as  {r  («)>0}.  For  a  given 

starting  point  within  the  feasible  region  the  iteration  process  is 

conducted  in  the  same  manner  as  in  the  unconstrained  problem.  However,  the 

iterate  is  now  deflected  away  from  the  constraint  boundaries  by  the  penalty  term 

* 

?(m).  When  the  process  has  converged  to  the  maximum  of  A  (within  the 

feasible  region)  the  process  is  restarted  from  the  point  of  convergence  using  a 

★ 

smaller  set  of  ic  values  in  the  objective  function  A  .  This  sequence  of 

.  *  * 
operations  is  repeated  until  there  is  no  appreciable  change  in  the  final  A 

values.  It  should  be  noted  that  only  those  constraints  in  danger  of  being 

violated  need  be  included  in  (4-10).  That  is,  we  might  set  a  number  of  the  k 

coefficients  to  zero  throughout  the  course  of  the  computations. 

In  the  present  problem  we  require  that  the  variance  covariance  matrices  r 

2 

and  T  be  positive  definite  and  a  >0.  As  r,  given  by  (2-17),  is  diagonal 

it  will  be  positive  definite  if  and  only  if  its  elements  are  positive.  And  T, 

given  by  (2-19),  will  be  positive  definite  if  and  only  if  the  diagonal  blocks 

©^,  i  =  1,2,... ,C  are  positive  definite.  It  is  convenient  to  specify  the 

constraints  on  the  {©^}  in  terms  of  their  discriminants.  The  mth 

discriminant  of  ©.,  denoted  as  6.  ,  is  defined  as  the  determinant  of  the 
l  lm 

upper  left  hand  submatrix  of  size  m.  (This  could  also  be  called  the  upper  left 


(4-9) 

(4-10) 
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hand  principal  minor  of  order  m.)  It  can  be  shown  (see  Hildebrand,  Reference 
28,  p.  51)  that  the  p  x  p  matrix  is  positive  definite  if  and  only  if 
>  0,  m  *  1,  2,...,p.  Consequently,  we  may  now  specify  the  required 
constraints  in  a  manner  that  is  consistent  with  the  r  notation  used  in  (4-10). 
These  are 


2 

a 


>  0, 


(4-11) 


Yk  >  0,  k  *  1,2,. ...K  (4-12) 

5im  >  0;  1  ”1»2»***»C*  m  *  1.2,. ..,p  .  (4-13) 

In  order  to  use  the  unconstrained  techniques  developed  earlier  in  this 

* 

section  for  the  purpose  of  maximizing  A  we  need  first  and  second  partial 

derivatives  of  t(w).  These  are  then  subtracted  from  the  corresponding 

derivatives  of  A  (or  X^)  to  form  the  derivatives  of  A*  (or  A^). 

Making  a  self-explanatory  change  in  the  k  notation  and  defining  c . „  as  the 

ltm 

cofactor  of  the  element  in  the  matrix  with  determinant  we  can 

write  the  derivatives  as 


(4-14) 

(4-15) 

^Hildebrand,  f.  B.,  Methods  of  Applied  Mathematics  (Englewood  Cliffs, 

NJ:  Prentice-Hall,  Inc.,  1952). 


i£_  -  -  k<o2)/o4 
3o2 


iJL_  ■  -  <(y  )/y2 
3YW  k  k 
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ie(5  )  c  /62 
im  i  tm  in 


(4-16)* 


JLS —  -  2  <(a2)/a6 

3a2 3a2 


(4-17) 


3  C  =  2  K(Y  )/y3 

3Yk3Yk  k  k 


(4-18) 


-Li _  =  2  V  <(5.  )  c2  /63 

3Tit3tit  m=l  in  itm  im 


(4-19) 


«  Z  «<«.  >  k  c  /«3  -  .8Ci..tn!  / 62  1 


(4-20) 


3  Yk  3<J2  3tj  ».  3a2  3Ylc3Y*  3Yk3nt  3Ths3Tit 


In  (4-21)  k  ^  i.  and  h  f  i.  For  a  value  of  p  equal  to  2  or  3  the  value  of 

3c^tm/3t£g  in  (4-20)  is  0  or  +1  respectively  with  the  signs  of  the 

latter  depending  on  the  detailed  descriptions  of  the  discriminant  functions 

define  a  cofactor  c.  as  zero  if  the  5.  discriminant  function  does  not 
itm  im 

contain 


♦Derivatives  of  a  determinant  j A j  may  be  obtained  from  its  expansion  in 
terms  of  cofactors  a|  *  2.  a^j  c^j,  where  cjj  is  the  cofactor  of 
ajj  in  A.  Hence,  3  A|/3ajj  *  c£j  (see  Searle,  Reference  20,  p.  86). 

20See  footnote  20  on  page  3-10. 
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4-3  OTHER  COMPUTATIONAL  CONSIDERATIONS 


Monitoring  the  progress  of  the  iteration  process  in  either  the 

,  ★ 
unconstrained  or  constrained  settings  will  require  the  evaluation  of  X  or 

X  .  For  this  purpose  we  write  them  as 

X*  =  -  N  log  (2iro2)  -  I  log  (|H|)  -  _L_  (£-Xu-Zv)  (4-22) 

2  2  2o2 

X#  *  -  ^N~Cp)  log  a2  -  I  log  (I  Hi  *1  x'h-1XI  )  -  _L_  y'(y-XTT-Zv),  (4-23) 

2  2  2o2  - 

which,  with  the  exception  of  the  determinants,  involve  quantities  that  are 
easily  computed  from  previous  results. 

We  can  obtain  the  determinants  |  H |  and  |X  H  X  |  in  forms  more  suitable 
for  computations  by  using  the  identities  (see  Searle,  Reference  20,  p.  96) 


pll 

Gl2l 

,  -1 

iGlll 

•  |G22  -  G21  G11  g12 1 

(4-24) 

LG21 

G22j 

pll 

Gl2] 

.  -1 

5 

lG22l 

1 G1 1  ~  g12  g22  g21 1 

,  (4-25) 

ig21 

G22j 

which  respectively  assume 

the  nonsingularity  of  and  G^- 

By  applying 

these  identities 

to  expand 

the 

determinant 

2GSee  footnote  20  on  page  3-10. 
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(D-1  +  Z'Z)  Z* 

Z  I 

we  find,  as  shown  by  Hemmerle  and  Hartley  (Reference  11),  that 

| H  |  S  |D-1+Z'Z|  •  j  D |  5  | B  j  •  j  D  |  .  (4-26) 

Now,  using  the  partitioned  form  of  B  given  in  (3-20)  and  the  above 
identities  we  can  write  | B  |  in  two  forms 

u'u+r-1  u'w 

|  B I  »  •  (  =  ju'u+r"1 1 »|  w'w+T“1-w'u(u'u+r~1)~1u,w  |  (4-27) 

w’u  w'w+T-1 

and 

u'u+r-1  u'w 

j  B  |  -  t  (  =  | w'w+T-1]  *|  u'u+r-1  -  u'w  (w'w+T-1  )-1w'u|  .  (4-28) 

w'u  w'w+T-1 

It  is  observed  that  the  right  most  determinants  in  these  equations  involve 

matrices  that  are  inverted  to  give  ¥22  and  in  (3-22a)  and  (3-24b). 

Since  the  determinant  of  a  matrix  can  in  many  cases  be  obtained  as  a  by-product 

of  the  inversion  algorithm,  e.g.,  a  Cholesky  decomposition  (Reference  11),  we 

can  easily  compute  |B|  by  either  (4-27)  or  (4-28),  depending  upon  which  set  of 

equations  (3-22a)  through  (3-24a)  or  (3-22b)  through  (3-24b)  are  used  to  invert 

B.  The  other  determinants  needed  to  compute  H  are  also  easily  obtained; 

|U  U+r  |  is  the  product  of  diagonal  elements  and  | D]  and  jW  W+T  |  are 

given  by  the  products  of  the  determinants  of  the  diagonal  blocks.  To  obtain 
*  "1  _ 

| X  H  X|  we  note  that  given  by  (3-15)  may  be  expressed  as 

^See  footnote  11  on  page  3-1. 
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AX1  -  |x'  [I-ZD(I+Z,ZD)_12'  ]x}_1  -  [x'h^X]"1,  (4-29) 

»  _i 

from  (3-7),  (3-18),  and  (3-39).  Hence,  ]X  H  X|  may  be  formed  as  a 
by-product  of  the  inversion  process  by  which  we  obtain  A^. 

As  a  final  note  concerning  the  numerical  methods  it  is  observed  that  many 
of  the  matrix  operations  required  to  invert  matrices  A  and  B  and  to  calculate 
the  derivatives  can  be  done  once  and  reused  in  the  subsequent  iterations.  For 
example  the  products  x'x,  x'u,  x'w,  u'u,  u'w,  w'w,  £*X,  Z>u»  Zw>  and  z'z 
involve  known  fixed  constants.  Furthermore,  since  the  various  equations  used  in 
the  calculations  depend  on  the  larger  (N  rows)  matrices  X,  U,  W,  and  £  through 
these  products  only,  storage  can  be  economized  by  computing  the  product  matrices 
in  advance  and  not  saving  X,  U,  W,  and  £.  It  is  also  pointed  out  that  in 
computing  the  trace  of  a  matrix  product,  an  operation  that  frequently  appears  in 
the  expressions  for  the  derivatives,  it  is  only  necessary  to  compute  the 
diagonal  elements  of  the  product.  For  this  reason  and  because  the  matrices 
appearing  in  these  traces  are  generally  of  modest  sizes,  the  computation  of  the 
derivatives  should  be  rapid  in  comparison  with  the  inversion  times  for  matrices 


A  and  B. 


CHAPTER  5 


APPLICATIONS  OF  THE  MODEL 


In  this  section  we  indicate  some  of  the  applications  for  which  the  model 
may  be  used.  No  attempt  is  made  to  be  exhaustive  in  this  effort;  rather,  we 
will  attempt  to  show  by  a  few  useful  examples  that  the  applications  are  wide 
ranging. 

5-1  DIRECT  RESULTS 


Results  obtained  directly  from  the  model  computations  are  the  ML  or  REML 

estimates  of  the  following  quantities  (again  we  show  just  the  ML  estimates): 

.  2 

the  general  experimental  error  variance  a  ;  the  gage  calibration  error 
variances  for  different  gage  classes  given  by 

=  Var(a*  )  m  y  ,  k  *  1,2,...,  K;  (5-1) 

k  km  k 

the  variances  and  covariances  of  the  performance  variations  of  different 
explosive  classes  given  by 


o* 


1,2,...,  C 


(5-2) 


the  similitude  equation  coefficients  for  different  explosive  classes  given  by 
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the  sample  values  of  the  gage  class  calibration  errors  a  and  explosive 
performance  variations  b  given  by 

(f 

and  the  error  vector  or  vector  of  residuals  given  by 

|=I-XG-Ua-w£  .  (5-5) 

In  these  expressions  and  below  the  submatrices  of  A  ^  are  computed  using  the 
ML  or  REML  estimates  of  the  variance  covariance  matrices. 

It  is  common  in  ordinary  regression  models  to  examine  the  model  residuals 
for  systematic  trends  that  would  indicate  inconsistencies  between  the  data  and 
the  model.  Various  techniques  in  this  regard  are  discussed  by  Draper  and  Smith 
(Reference  29)  and  by  Seber  (Reference  7).  In  the  present  mixed  model  such 
examinations  can  be  extended  to  include  the  estimates  of  all  the  random 
effects.  Tests  of  the  distributional  assumptions  can  be  made  directly  on  the 
values  of  £,  £,  and  £  using,  for  example,  the  methods  summarized  in  Mehrotra  and 
Michalek  (Reference  30).  The  identification  of  faulty  gages  or  of  anomalously 
performing  explosive  charges  should  be  immediately  evident  upon  examination  of  £ 

A 

and  b. 


[A*  X*  +  A  z']£  ;  (5-4) 

12  22 


29Draper,  N.  R. ,  and  Smith,  H. ,  Applied  Regression  Analysis  (New  York:  John 
Wiley  &  Sons,  1966). 

?See  footnote  7  on  page  2-6. 

■^Mehrotra,  k.  G.,  and  Michalek,  J.  E.,  "Tests  for  Univariate  and  Multivariate 
Normality,"  RADC-TR-76-140 ,  May  1976. 
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5-2  VARIANCE  COVARIANCE  MATRIX  FOR  u,  a,  AND  £ 


From  (3-12)  u,  and  v 


(a ' ,b  ' ) ' 


may  be  written  as 


where 


CRii 

■  H3  ■ 

A]^X  +A^2Z 

—i  i  —  i 

LR2J 

A  X  +A  Z 

L  12  22  . 

(5-7) 


Here  the  components  of  A  *  are  evaluated  from  the  ML  or  REML  estimates  of  the 

2  -1  itt 

variance  covariance  component  ratios  6.  and  a  .  Since  A  and 
are  complicated  functions  of  £  there  is  no  known  exact  expression  for  the 
variance  covariance  matrix  of  jl,  a,  and  {?.  Nevertheless,  it  is  a  common 
practice  in  this  situation  (e.g.,  see  References  14,  p.  205  and  15,  p.  35)  to 
obtain  an  approximate  variance  covariance  matrix  by  treating  A  1  as  though  it 
were  calculated  from  the  true  fixed  values  of  6.  Equivalently,  this  presumes 
knowledge  of.  the  true  value  of  H.  Equation  (5-6)  then  becomes  a  linear  function 
of  £  and  the  variance  covariance  matrix  is  easily  written  down  as  (see  footnote 
on  page  2-6). 


Var 


Rll 

(I+ZDZ' ) 

r2J 


a 


2 


Ri (I+ZDZ ' )Ri 
R2( I+ZDZ ' )R  1 


Rj ( I+ZDZ ' )R 2 
R2( I+ZDZ'  >112 


(5-8) 


^See  footnote  14  on  page  3-2. 
l^See  footnote  15  on  page  3-2. 
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Following  Henderson  (Reference'  31,  Appendix  A)  this  can  be  simplified  by  use  of 
the  relation 


rv 

[x,z]  * 

ah 

A12 

(rx’x 

x'z 

— « 

1  1 

Lr2J 

A 

L 12 

A 

22_ 

(Lz  X 

D-1+Z  Z. 

0 

0 


~1  O' 

0  A  D-1 

m 

12 

.0  I 

0  A  D-1 

22  J 

Hence,  R^X  *  I 

RlZ  “  -  Ai2 
R2X  »  0 

R2Z  =  I  -  A22D~^  . 


(5-9) 


(5-10) 


By  substituting  (5-7)  and  (5-10)  into  (5-8)  one  -readily  obtains 


A11 

0 


0 

D-A22 


a 


2  . 


(5-11) 


An  estimate  of  (5-11)  can  be  obtained  by  substituting  the  ML  or  REML  estimates 
of  D,  A^,  A^2»  and  cr  .  The  result  Var(j£)  *  A^a  could  have  also  been 
derived  as  done  by  Corbeil  and  Searle  (Reference  15)  from  (3-5)  and  (4-29). 


In  Reference  14  (p.  205)  Harville  has  pointed  out  that  (5-11),  because  H  is 
assumed  to  be  known,  will  tend  to  underpredict  the  dispersions  of  u  and 


^Henderson,  C.  R. ,  "Best  Liner  Unbiased  Estimation  and  Prediction  Under  a 
Selection  Model,"  Biometrics,  Vol.  31,  p.  423. 


15see 

footnote 

15 

on  page 

3-2 

14See 

footnote 

14 

on  page 
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3 


I 


i 


H 

►  * 


s 


I 


►** 

>  v 

k 

r*- 

►  . 


ri 

»•  - 


i\ 

[» 


A  i  Ai  i 

(a  ,b  )  but  suggests  also  that  "the  downward  bias  may,  at  least  in  some 
instances,  be  so  small  as  to  be  unimportant."  Ouantif ication  of  the  accuracy  of 
this  result  (as  well  as  of  other  results  derived  below  in  a  similar  manner)  is 
possible  by  means  of  Monte  Carlo  simulation  studies.  Such  studies  will  be 
pursued  upon  completion  of  the  model  coding  and  published  in  a  subsequent  report. 


5-3  LARGE  SAMPLE  VARIANCE  COVARIANCE  MATRICES  FOR  ESTIMATORS  OF  COMPONENTS  OF 
VARIANCE  AND  COMPONENTS  OF  VARIANCE  RATIOS 


A  well  known  result  from  the  theory  of  maximum  likelihood  (e.g.,  see 
Kendall  and  Stuart,  Reference  32.)  is  that  the  asymptotic  (N-»-«>)  variance 
covariance  matrix  of  the  vector  of  parameter  estimates  is  given  by  the  inverse 
of  the  information  matrix  J,  which  is  the  negative  of  the  matrix  of  expected 
values  of  the  second  derivatives  of  the  log-likelihood.  In  our  case  with  (4-1) 
as  the  vector  of  parameters  we  have 

J=-E<X),  (5-12) 


where  J6  is  given  by  (4-4).  Hence,  denoting  the  vector  of  ML  estimates  as 


A  /  ^9  A  A  A  /N  \  1 

j£  -  '  >  Yl  »  •  •  •  >  TR*  Til  »  *  *  *  T1  ic a  •  •  •  ,  Tq].  ,  .  .  •  ,  T£  v) 


(5-13) 


we  can  write 

Var(w)  >  J*1, 


(5-14) 


^Kendall,  M.  G. ,  and  Stuart,  A.,  The  Advanced  Theory  of  Statistics  (London: 
Charles  Griffin  &  Co.,  Ltd.,  1973),  Vol.  2. 
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where  equality  is  achieved  in  the  limit  as  N-*-®.  This  lower  bound  is,  however, 
often  used  as  an  approximation  for  Var(u).  An  estimate  of  Var(u)  may  thus 
be  obtained  from  the  final  iterated  value  of  ECJ0. 


Often  it  may  be  of  more  interest  to  know  the  variance  covariance  matrix  of 
the  estimators  of  the  components  of  variance  rather  than  of  the  components  of 
variance  ratios.  Suppose  this  is  denoted  as 


^  /  /S0  *9  <*9  ^  A  /V  /A  v  1 

^  -  \  O  j  •  •  •  |  0  jU  j  •  •  •  V)  j  •  •  *  j  \J  j  •  •  *  j  l)  / 

~  1  K  11  1  IT  Cl  Cir 


(5-15) 


where  as  in  (5-1)  a2  =  y  a2,  k  *  1,...,K  and  £  =  £  o2,i  =  1,...,C,  j  = 

k  k  ij  ij 

l,...,if.  Corbeil  and  Searle  (Reference  15)  have  shown  that  the  relationship 
between  Var(<?)  and  Var(w)  is  given  by 


Var  (£) 


0  Var(w) 


n 


(5-16) 


where  n  is  the  Jacobian  matrix  for  the  transformation  from  to  w.  This 
is  easily  shown  to  be 


(5-17) 


2  '  ' 

where  m-  is  defined  by  =  (a  ,w_)  .  Hence,  the  large  sample  value  of 

Var(£)  can  be  obtained  by  substituting  J  *  for  Var(5)  in  (5-16). 


l^See  footnote  15  on  page  3-2. 
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5-4  PREDICTION  OF  THE  PERFORMANCE  IN  A  PAST  TEST 


In  our  discussion  of  Che  basic  model  for  a  single  observation,  Section  2-1, 

we  let  f^j(x)  denote  the  value  of  the  principal  quantity  of  interest  in  the 

shot  of  the  ith  explosive  class  at  some  reduced  travel  distance  x  from  the 

charge.  We  assume  f„(x)  is  any  1  to  1  transformation,  such  as  the  logarithm, 

of  a  possibly  reduced  performance  index  or  effectiveness  factor  that  has 

physical  meaning  and  significance.  In  the  light  of  past  discussions  we  will  now 

use  a  somewhat  more  explicit  notation  and  let  f..(x)  denote  the  realization  of 

ij 

ic 

the  random  variable  f„(x)  that  in  accord  with  (2-7)  may  be  expressed  as 


=  i'jii  +  ±'£ij 


(5-18) 


Here  denotes  some  vector  function  of  the  arbitrary  distance  x  analogous  to 
■^ijkm  ^escri^ec*  earlier.  Following  suit,  we  write  the  unknown  realization 

"fa 

of  f.j(x),  that  is  the  sought  after  quantity,  as 


fij(*')  =  1  M.£  +  IJij  , 


and  the  ML  estimator  of  f..(x)  as 

ij 


Fij(x)  =  £  £.  +  •  g-j 


(5-20) 


A  . 

Her?  3-  •  and  3.  ■  are  the  subvectors  of  the  vectors  o  and  b  that  correspond 
-ij  ^ij  y 

jth  shot  of  the  ith  explosive  class. 
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Confidence  limits  for  f^j(x)  may  be  obtained  under  the  assumption  of  a 
known  variance  covariance  ratios  matrix  H.  The  limits  so  obtained  may  be 
reasonably  accurate  when  sample  sizes  are  large,  but  when  samples  are  only 
modestly  sized  the  confidence  interval  will  correspond  to  a  confidence 
coefficient  that  is  smaller  than  the  one  specified.  Nevertheless,  it  may  be 
possible  to  attach  a  more  realistic  confidence  coefficient  to  the  interval  by  a 
Monte  Carlo  simulation  procedure.  Hence,  we  include  the  theory  in  this  report. 


Dropping  the  argument  indicating  the  explicit  dependency  on  x,  we  can 
define  a  vector  y  such  that 


-ij 


(5-21) 


Comparing  this  with  (5-19)  it  is  obvious  that  consists  of  a  column  of  zeros 
imbeded  with  two  vectors  located  in  such  a  way  as  to  extract  the  and 
0.  .  subvectors  from  and  v.  Similarly  we  can  write 


and 


(5-22) 


(5-23) 


The  development  of  confidence  limits  for  f„  can  proceed  in  a  manner 
similar  to  that  found  in  a  more  general  study  due  to  Farville  (Reference  33). 
Under  the  assumption  of  known  variance  covariance  ratios,  y  is  unbiased  and  it 


33Harville,  D.  A.,  "Confidence  Intervals  and  Sets  for  Linear  Combination  of 
Fixed  and  Random  Effects,"  Biometrics ,  Vol.  32,  1976,  p.  403. 
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can  be  shown  (see  Henderson,  Reference  31  for  an  explicit  derivation)  that  the 


random  variable  f . . -f . .  is  normally  distributed  with  a  zero  mean  and 

ij  ij 


1  -1  2 

variance  f  A  Ta  ;  hence 


(  f i j-f i j )/(l'A"1fo2)1/2  ~  N(0 , 1 ) 


(5- 


Furthermore , 

'a2  =  Z'e/(N-Cp)  (5- 

2  2  2 
is  an  unbiased  estimator  of  a  and  it  can  be  shown  that  (N-CP)or  /a 

2 

is  x  distributed  with  (N-Cp)  degrees  of  freedom  and  is  independent  of 
(5-23).  From  these  results  it  follows  that 


f.  f*. 

_LJ _ iJ _  ~  t  (5- 

(^A~lpl/2^  N-Cp 

i.e.,  has  a  t  distribution  with  N-Cp  degrees  of  freedom.  Then,  from  the 
symmetry  of  the  t  distribution,  we  have 

Pr[  |Cij-f!j  |  <  (xVl?)l/2  ?  t6/2jN_Cp]  -  1  -  «  ,  (5- 

where  t(5/2  N-Cp  t^ie  ^06/2  percentile  of  the  tN_Cp  distribution. 

Hence,  for  a  particular  realization  we  find  that 

fij  *  (i'a-1!)172  7  t6/2jN_Cp  (5- 


31 


See  footnote  31  on  page  5-4 
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are  1-5  confidence  limits  for  f^j  (termed  "unconditional"  confidence  limits 
by  Harville  in  Reference  33).  These,  of  course,  are  valid  for  any  values  of  the 
regressor  variable  x. 


1  -i 

For  computation  purposes  it  is  useful  to  express  !  A  f  m  the 
following  manner 

A~1  +  2*_  A^2^Jli>^ij )  +  422^ij^  £.  • 

Here  and  c*enote  P  x  ?  submatrices  of  A  having 

rows  and  columns  corresponding  to  and  0.  ^  respectively  and  the 

submatrix  A,  ~(u* >#• • )  corresponds  to  the  u-  rows  and  the  8.. 

12  — l  — 1  j  r  —i  — 1J 

columns . 


(5-29) 


If  a  logarithmic  transformation  of  the  response  variable  has  been  employed 
it  may  be  of  interest  to  estimate  and  obtain  confidence  limits  for  the  antilog 
of  fjj.  Suppose  the  common  logarithm  was  used.  Then  the  ML  estimator  of  the 
(possibly  reduced)  performance  index  in  the  jth  shot  of  the  ith 

explosive  class  would  be 


lij  "  10 


(5-3C 


and  from  (5-27)  confidence  limits  for  Tjj  would  be  given  by 
[fij  ±  (t'a-W)1/2  a  t«/2,N-Cp] 


(5-31 


33see  footnote  33  on  page  5-8. 
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5-5  PREDICTION  OF  THE  PERFORMANCE  IN  A  FUTURE  TEST 


Continuing  the  notation  of  the  previous  section,  we  now  consider  the 
prediction  of  the  performance  index  in  a  future  test  at  an  arbitrary  distance  x 
from  the  charge.  We  assume  the  charge  is  of  the  ith  explosive  class  where 
1  £  i  _<  C.  To  avoid  confusion  with  past  quantities  we  will  indicate  the 
future  value  by  f.*,  i.e.,  an  asterick  in  place  of  the  subscript  denoting  a 
particular  shot. 


In  analogy  with  (5-18)  we  let  f.*  be  the  realization  of  the  random 

1* 


variable 


f£*  *  i  ui  +  2.&X*  • 


(5-32) 


where  0.*  is  independent  of  past  quantities  and  therefore  of  v  and 

Under  the  assumption  of  a  known  H,  an  unbiased  estimator  of  f.  is  given  by 

1* 


f.  .  *  9  u.  . 

i*  —  —i 


(5-33) 


To  obtain  confidence  limits  for  f ^ under  the  same  assumption  consider 

*  it 

the  random  variable  f. .-  f. ..  We  have 

l*  l* 


E  (f,*-f.-*)  =  *  u-  -  #  u-  =  0 

i7'  —  — l  —  — i 


Var(?i*-f**)  =  ♦  (Au(i£)  +  e.  )  to2  , 


so  that  f.-*-f.-*~N(0,  *  (A, ,  ( p • )  +0.)  9a2). 


(5-34) 


(5-35) 


(5-36) 


5-11 


we  have 


'fg 

Now,  since  a  of  (5-24)  is  also  independent  of  f. 

1" 

1  -f* 

**  _  _t  (5-37) 

[  $_'  (An (  Mi )  +  ei)*]1/2a  N‘CP 

This  allows  us  to  make  the  probability  statement 

Pr{|fi*  -  ff*  |  <  [*,(A11(J!i)  +  a  t{/2>N_Cpl  -  1  -  «  •  (5-38) 

Hence,  1-6  confidence  limits  for  a  particular  future  realization  f^* 
ard  given  by 

fi*  -  t  (All(lii)  +  o  t  6/2 ,N-Cp  *  (5-39) 

As  emphasized  in  the  previous  section,  these  limits'  will  tend  to  underpredict 
the  size  of  the  confidence  interval  because  we  have  assumed  the  estimated  value 
of  H  to  be  the  true  value. 

If  a  common  logarithmic  transformation  of  the  data  has  been  used,  the  ML 
f. 

estimate  of  T^*  =  10  x*  can  be  obtained  from 

A 

*  f  • 

T£*  *  10  ,  (5-40) 

and  approximate  confidence  limits  for  T£*  from 

(f£*  ±  [♦'(Ai1(jj£)  +  ) # ] 1  / 2  ^  t5/2,N-Cp} 

10 


(5-41) 
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